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Abstract 

This paper gives an in-depth study of a multiple-antenna wireless communication 
scenario in which a weak signal received at an intermediate relay station is am- 
plified and then forwarded to the final destination. The key quantity determining 
system performance is the statistical properties of the signal-to-noise ratio (SNR) 
y at the destination. Under certain assumptions on the encoding structure, recent 
work has characterized the SNR distribution through its moment generating func- 
tion, in terms of a certain Hankel determinant generated via a deformed Laguerre 
weight. Here, we employ two different methods to describe the Hankel determi- 
nant. First, we make use of ladder operators satisfied by orthogonal polynomials 
to give an exact characterization in terms of a "double-time" Painleve differential 
equation, which reduces to Painleve V under certain limits. Second, we employ 
Dyson's Coulomb Fluid method to derive a closed form approximation for the 
Hankel determinant. The two characterizations are used to derive closed-form 
expressions for the cumulants of y, and to compute performance quantities of 
engineering interest. 
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1. Introduction 

Over the past decade, multiple-input multiple-output (MIMO) antenna sys- 
tems employing space-time coding have revolutionized the wireless industry. Such 
systems are well-known to offer substantial benefits in terms of channel capacity, 
as well as improved diversity and link reliability, and as such these techniques are 
being incorporated into a range of emerging industry standards. More recently, 
relaying strategies have been also proposed as a means of further improving the 
performance of space-time coded MIMO networks. Such methods are particularly 
effective for improving the data transmission quality of users who are located near 
the periphery of a communication cell. 

Various MIMO relaying strategies have been proposed in the literature; see, 
e.g., [1-6], offering different trade-offs in various factors such as performance, 
complexity, channel estimation requirements, and feedback requirements. In [7], 
a communication technique was considered which employed a form of orthog- 
onal space-time block coding (OSTBC), along with non-coherent amplify-and- 
forward (AF) processing at the multi-antenna relay. This method has the advan- 
tage of operating with only low complexity, requiring only linear processing at 
all terminals, achieving high diversity, and not requiring any short-term channel 
information at either the relays or the transmitter. For this system, an important 
performance measure — the so-called outage probability — which gives a funda- 
mental probabilistic performance measure over random communication channels, 
was considered in [7] and [8], based on deriving an expression for the moment 
generating function of the received signal to noise ratio (SNR). The outage prob- 
ability could then be computed by performing a subsequent numerical Laplace 
Transform inversion. The exact moment generating function results presented in 
[7] and [8], however, are quite complicated functions involving a particular Han- 
kel determinant, and they yield little insight. Moreover, for all but small system 
configurations, the expressions are somewhat difficult to compute, particularly 
when implementing the Laplace inversion. 

In this paper, we significantly expand and elaborate upon the results of [7] 
and [8] by employing analytical tools from random matrix theory and statistical 
physics. Technically, the challenge is to appropriately characterize the Hankel 
determinant which arises in the moment generating function computation. In gen- 
eral, the Hankel determinant takes the form, 




(1.1) 



2 



generated from the moments of a certain weight function w(x), < x < oo, 

oo 

fit := J" x k w(x)dx , k = 0, 1, 2, 

o 

Note that the above is merely one of several possible representations for D n , with 
equivalent forms regularly found in the fields of statistical physics and random 
matrix theory. 

For the problem at hand, we are faced with a Hankel determinant that is gen- 
erated from a weight function of the form, 

w AF (x, T, t) = x a e~ x (^-t , < x < oo, (1.2) 
\T + xl 

which is a "two-time" deformation of the classical Laguerre weight, x a e~ x . The 
parameters in this weight, satisfy the following conditions: 

a>-\, T := — — , t>0, c>0, N s > 0, < s < oo. (1.3) 
1 + cs 

For our problem both a and N s are integers, where a := \Nr - Nd\ is the abso- 
lute difference between the number of relay and destination antennas, Nr and No 
respectively 1 . 

We will apply two different methods to characterize this Hankel determinant. 
For the first method, we will derive new exact expressions for D n [w af] by employ- 
ing the theory of orthogonal polynomials associated with the weight waf(*> T, t) 
and their corresponding ladder operators. For such methods, extensive literature 
exists (see, e.g., [10-20] for their use in applications involving unitary matrix 
ensembles), though in the context of information theory and communications the 
techniques have only very recently been introduced by the first and third authors in 
[13]. For the second method, we derive an approximation for the Hankel determi- 
nant using the general linear statistics theorem obtained in [21], based on Dyson's 
Coulomb Fluid models [15, 22-25]. These results are essentially the Hankel ana- 
logue of asymptotic results for Toeplitz determinants [26]. Once again, these 
techniques have been used extensively, particularly in statistical physics, though 
only very recently have they been applied to address problems in communications 
and information theory [13, 27]. 



'The parameters a and N s are determined by the model, but in fact can be extended to take 
non-integer values such that a > -1 and N s > [9]. 
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It is important to note that in addition to the two methodologies advocated 
above, there exists other integrable systems approaches which can be used for 
characterizing Hankel determinants. For example, D n may be written in an equiv- 
alent matrix integral formulation, from which a 'deform-and-study' or isomon- 
odromic deformation approach may be adopted. Essentially, this idea involves 
embedding D n into a more general theory of the r-function [28, 29], using bilin- 
ear identities and linear Virasoro constraints. For details, see [30-32]. Yet another 
approach is to characterize D n through the use of Fredholm determinants as em- 
ployed by Tracy and Widom [10]. Both these exact, non-perturbative integrable 
systems methods have their own specific advantages and disadvantages, and in 
most cases lead to non-linear ordinary differential or partial differential equations 
(ODE/PDEs) satisfied by the Hankel determinant. However, these equations are 
usually of higher order 2 , from which first integrals have to be found to reduce to 
a second order ODE. The advantage of the ladder operator method is that closed- 
form second order equations are directly obtained for a quantity related to the 
Hankel determinant, bypassing the need to find a first integral, or evaluate any 
multiple integral. 

Based on the derived exact and asymptotic (large n) approximations for the 
Hankel determinant D„, and thus corresponding characterizations for the moment 
generating function of interest, we investigate the distribution of the received SNR 
of OSTBC MEMO systems with AF relaying, and also the error performance. In 
particular, the Coulomb Fluid representation for the moment generating function 
of the received SNR is shown to yield extremely accurate approximations for the 
error performance, even when the system dimensions are particularly small. We 
also employ our analytical results to compute closed-form expressions for the 
cumulants of the received SNR, first via the Coulomb Fluid method, and then pre- 
senting a refined analysis based on Painleve equations. Subsequently, we give an 
asymptotic characterization of the moment generating function, valid for scenarios 
for which the average received SNR is high, deriving key quantities of interest to 
communication engineers, including the so-called diversity order and array gain. 
These results are, once again, established via the Coulomb Fluid approximation, 
and subsequently validated with the help of a Painleve characterization. 



2 For Painleve systems, equations of Chazy type usually appear [33]. 
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1.1. Amplify and Forward Wireless Relay Model 

Here we briefly recall the background for the model, having been developed 
in [7]. The dual hop MIMO communication system features a source, relay and 
destination terminal, having N s , Nr and No antennas respectively. In the process 
of transmitting a signal, each transmission period is divided into two time slots. 
In the first time slot, the source transmits to the relay. The relay then amplifies 
its received signal subject to an average power constraint, prior to transmitting 
the amplified signal to the destination terminal during the second time-slot. We 
assume that the source and destination terminals are sufficiently separated such 
that the direct link between them is negligible (i.e., all communication is done via 
the relay). 

Let Hi e C NrXNs and H2 e C n ° xNr represent the channel matrices between 
the source and relay, and the relay and destination terminals respectively. Each 
channel matrix is assumed to have uncorrelated elements distributed as 3 CN(0, 1). 
The destination is assumed to have perfect knowledge of H 2 and either Hi or the 
cascaded channel H 2 Hi, while the relay and source terminals have no knowledge 
of these. 

A situation is considered where the source terminal transmits data using a 
method called OSTBC encoding [34]. In this situation, groups of independent and 
identically distributed complex Gaussian random variables (referred to as infor- 
mation "symbols") Sj, i = 1, . . .,N are assigned via a special codeword mapping 
to a row orthogonal matrix X = (xi, . . . , Xn p ) e C NsXNp satisfying the power con- 
straint E [||x t || 2 ] = y, where N P is the number of symbol periods used to send each 
codeword. 

Since it takes Np symbol periods to transmit N symbols, the coding rate is then 
defined as 

N 

R = — . (1.4) 

N P V ' 

The received signal matrix at the relay terminal at the end of the first time slot, 
Y = (yi, . . . , y Np ) G C^ x ^, is given by 

Y = HiX + N, (1.5) 



(T 



3 The notation CN(/u, o 2 ) represents a complex Gaussian distribution with mean fi and variance 

2 
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where N e C NrXNp has uncorrelated entries distributed as CN(0, 1), representing 
normalized noise samples at the relay terminal. The relay then amplifies the signal 
it has received by a constant gain matrix G = fll^, where 

a 2 = - . (1.6) 

a+7)N R 

In the above, b is the total power constraint imposed at the relay, i.e., E [jIGy^ll 2 ] < 
b, whilst y represents the average received SNR at the relay. The received signal 
R e C Nc,xNp at the destination terminal at the end of the second time slot is then 
given by 

R = aH 2 Y + W = 5H 2 HiX + «H 2 N + W, (1.7) 

where W e C Nc,xNp has uncorrelated entries distributed as CN(0, 1), representing 
normalized noise samples at the destination terminal. 

Next, the receiver applies the linear (noise whitening) operation, 

R = K 1 /2 R, K := a 2 H 2 H^ + I Nd , (1.8) 

to yield the equivalent input-output model 

R = HX + N, (1.9) 

where H = <3K~ 1/2 H 2 H 1 and N e C NoXNp has uncorrelated entries distributed as 
CN(0, 1). 

Based on the above relationship, standard linear OSTBC decoding can be ap- 
plied (see, e.g., [34]). This results in decomposing the matrix model (1.9) into a 
set of parallel non-interacting single-input single-output relationships given by 

si = ||H|| F 5, + T] U l=l,...,N, 

where rji is distributed as CN (0, 1), with ||H|| f representing the Frobenius norm 
(or matrix norm) of H. 

The quantity that we are interested in, the instantaneous SNR for the Zth sym- 
bol yi, can then be written as 

71 = n fi ii^M 2 ] = ^.(i%)^ Tr ( H ^ 2K " lH2Hi )- (L10) 

Since the right-hand side is independent of /, we may drop the / subscript and 
denote the instantaneous SNR as y without loss of generality. 
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1.2. Wireless Communication Performance Measures 

Here we recall some basic measures which are employed to quantify the per- 
formance of wireless communication systems. One of the most common mea- 
sures is the so-called symbol error rate (SER), which quantifies the rate in which 
the transmitted symbols (or signals) are detected incorrectly at the receiver. 

For signals which are designed using standard M-ary phase shift keying (MPSK) 
digital modulation formats, the phase of a transmitted signal is varied to convey 
information, where M e {2, 4, 8, 16, . . .} represents the number of possible signal 
phases. The SER can be expressed as follows [35] 

Pmpsk = - f aJ^FW (1.11) 
n J \ sin 2 6 j 
o 

where Aty(-) is the moment generating function of the instantaneous SNR y at the 
receiver, and © = n(M - l)/M and gMPSK = sin 2 (7r/M) are modulation-specific 
constants. As also presented in [36], the SER may be approximated in terms of 
Aty(-) but without the integral, via the following expression: 

(1.12) 

In addition to the SER, another useful quantity is the so-called amount of 
fading (AoF), which serves to quantify the degree of fading (i.e., the level of 
randomness) in the wireless channel, and is expressed directly in terms of the 
cumulants of y. Specifically, the AoF is defined as follows: 

AoF = k 2 /k 2 , (1.13) 

where K\ and k 2 are the mean and variance of y respectively. As shown in [37], 
for example, this quantity can be used to describe the achievable capacity of the 
wireless link when the average SNR is low. 

For our AF relaying system under consideration, it is clear that a major chal- 
lenge is to characterize the moment generating function and also the cumulants of 
the instantaneous SNR given in (1.10). This is the key focus of the paper. 

1.3. Statistical Characterization of the SNR y 

Here we introduce the moment generating function and cumulants of interest, 
and pose the key mathematical problems to be dealt with in the remainder of the 
paper. 
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Defining the positive integers, 

m := max(N R ,N D ), n = min(N R ,N D ), 

where n + due to physical considerations, it was shown in [7] that the mo- 
ment generating function of the instantaneous SNR (1.10) may be expressed as 
the multiple integral, 



h f n (*j - xif ft *r n e-» ( » ) s d Xk 

" miy \<i<j<n k=l \ l+a ( 1 + M7^/ 

M y (s) = —-- . (1.14) 

h J n (x 7 - - x ; ) 2 n 

[0,oo )» 1<'<7<'1 k=l 



From the above, we may then compute the 1 th cumulant of y via 



17' 



(1.15) 



,v=0 



Let 



a := m - n, (1-16) 

t:=zs, (1.17) 
a 2 

c:=i (1.18) 

r = r(,): =(i7S)- (U9) 

We now consider the moment generating function (1.14) as a function of two 
variables, (T, t), or (5,0. since T = t/(l + cs). We may then write the moment 
generating function (1.14) as 

N i / n - ^ ft (£|f 

M r (r,r)= - . (1.20) 

h I n (x 7 -xo 2 n^-^ 

[0,oo)" !<!<;<« k=i 
We see that this involves the weight with the parameters T and t 

w AF (x,T,t) = jfe-'^Y^' , (1.21) 
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which is a deformation of the classical generalized Laguerre weight 

w^l g (x) = x a e~ x , < x < oo. (1.22) 

Note that w AF (x, t, t) = w ( ^ g (x). 

The multiple integral representation (1.20) is expressed as a ratio of Hankel 
determinants: 

/TV* D w [ WAP (-,r,Q] 
AWr '° U D n [<(.)] ' (L23) 

T nN s det( nt+j(T,fj). . 

■7 ~7 F 2 - (1 ' 24) 

where ///(7\ r) is the /th moment of the weight 

oo 

fj.j(T,t)= Jx J w AF (x,T,t) dx, j = 0,1,2,... 
o 

For A^. e N, the moments of the weight w AF are expressed in terms of the Kummer 
function of the second kind, 

fXj(T, t) = t a+ i +l T(a + j + 1) ^ P\ U(a + j + 1, a + j + 2 - k, T). 

(1.25) 

The Hankel determinant for T = t, namely, 

D n [ WAF (;t,t)] = D n [w^(-)l 

can be regarded as a normalization constant, so that M y (t, = 1, and its closed 
form expression is well-known (see [38]). 

n-\ 

D n [w AF (-, t, 0] = Y\ T(a + i + l)T(i + 1), 

_ G(n + l)G(rc + a + 1) 
G(a + 1) 
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where G(z) is the Barnes G-function defined by G(z+ 1) = Y(z)G(z) with G(l) = 1 . 

From a simple change of variables, the cumulants can be calculated by trans- 
forming (1.15) using (1.19) as 



T 2 d " 



K ' = 'Tjf' x °z M ^ T ^ 



(1-27) 

T=t 



Note that an equivalent expression to the result given in (1.25) has been pre- 
viously reported in [7]. This result, whilst having numerical computation advan- 
tages when the system dimensions are small, has a number of shortcomings. Most 
importantly, its complexity does not lead to useful insight into the characteristics 
of the probability distribution of the SNR y, and it does not clearly reveal the de- 
pendence on the system parameters. Moreover, in this form, the expression is not 
amenable to asymptotic analysis, e.g., as the number of antennas grow sufficiently 
large, and in such cases its numerical computation becomes unwieldy. 

In the following, we will seek alternative simplified representations with the 
aim of overcoming these shortcomings. The key challenge is to characterize the 
Hankel determinant in (1.23). 

Remark 1. It is of interest to mention here that for the special case a = 0, the 
Hankel determinant (1.23) is directly related to the shot-noise moment generating 
function of a disordered quantum conductor. This was investigated in [39] using 
the Coulomb Fluid method. Specifically, the two moment generating functions are 
related by 

M y (T,t)\ a=Q = M ^j^ T ^\ ( i.28) 
where z '■= t/T and Mshot-Noise(T, z) is the Hankel determinant generated 

w(x,T,z) = e- Tx (P^) Ns , 0<x<oo. (1.29) 



1 + x 

Of course, T, z and N s have different interpretations in this situation. 

We mention here that the operator theory approach of [40] justifies the Coulomb 
Fluid results obtained in [41 ] on the shot-noise problem. 

1.4. Alternative Characterization of the SNR y 

An alternative characterization for the moment generating function that will 
also prove to be useful is derived as follows. From the definitions of m and n, 
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there are two possible choices for the parameter t, 



1 (1+7)Nr 



' (l+y)« 
b 

( \+y)nr 



where 



r := m/n, 



N R < N D , 
N R > N D . 

(1.30) 



We first consider the sub-case N R < N D . To this end, starting with (1.20) and 
with the change of variables, x t — » nx h i = 1, . . . , n, we obtain 4 



My(T',t') 



where 



( \nN s 



Jj / n dx, exp 2 2 log \xj - x k \ - n £ (jt; - /? log *,) - , N s log ) 

[0 «>)'' '=1 L l<j<k<n j=\ \ J / 



ttt J n^/exp 2 2 log |*y - J*| - n £ ( - /3 log JC7I 
'[0,oo)» 1=1 I l<j<k<n 7=1 



Note that in terms of 5, T can be written as 



(1.31) 



t' 


t 

n 


(1.32) 


r 


T 
n 


(1.33) 


P 


m 
n 


(1.34) 



V 



Equivalently, we can write (1.31) as 



where 



Z n (T',t') := D n (nT',nt') 



1 + cs 



M 7 (T',f) = \- 



t \ nN s 



z n (.r,n 



f Z n (t',f) 



(1.35) 



(1.36) 



^ J exp -0(^,...,^)-|]^log|^-^J f]^-, (1.37) 

[0,co)» 1=1 1 -) =1 



For the sake of brevity, instead of defining a new function, we write M y (T' ', f ) in place of 
My(nT\nt'). 
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and 

n 

0(x!,...,x„) := -2 ^ log|x^-x fc | + n^(x / -j01ogx / ) . (1.38) 

1 < j<k<n i= 1 

Remark 2. We introduce the variables T and t' in order to account for the n- 
dependence of the variables T and t. This is important since the above representa- 
tion will be useful for deriving a large n approximation for the moment generating 
function based on the Coulomb Fluid linear statistics approach in Section 3. 

For the sub-case N R > N D , we have from (1.30) that t is instead given by 



b 



(1.39) 



where fi = - - 1. Hence equations (1.36)— (1 .38) and the results of Sections 3-5 
are valid for the case N R > N D upon transforming t' and T' to and 7^ via 

t' =: (l+jS)^ (1.40) 

and 

r =: (\+fi)T f (1.41) 
respectively. In this case, we may write in terms of the variable s as 

r f = . (1.42) 

1 + cs 

2. Painleve Characterization via the Ladder Operator Framework 

Using the ladder operator framework, and treating T and t as independent 
variables, the result below is proved in Appendix A. This gives a PDE satisfied 
by log M 7 (T, t) in the variables T and t. 

Theorem 1. Let the quantity H n (T, t) be defined through the Hankel determinant 
D n (T, t) as 

H n (T,t) := (Td T + td t )logD n (T,t) (2.1) 
= (Td T + td t ) log My(T,t), (2.2) 
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where the second equality follows from (1.23). Then H„(T, t) satisfies the follow- 
ing PDE 5 : 

H n = -2(d T H n )(d t H n ) + {In -N s + a + T)d T H n + (2n + N s + a + t)d t H n 

\(Td 2 TT + td 2 Tt )H n ]\(Td 2 Tt + tdl)H n ]+MH n )A 2 (H n ) 

+A l (H n )+A 2 (H n ) - V - ± , 

2[T(d T H n ) + t(d t H n ) -H n + n(n + a)] 

(2.3) 

where 

{Ai(H n )} 2 = ((Td 2 TT + td 2 Tt )H n f 

+4(T(d T H n ) + t(d t H n ) -H n + n{n + a))(d T H n )(d T H n - N s ), (2.4) 

and 

{A 2 (H n )\ 2 - ((Td 2 t + td 2 t )H n f 

+4(T(d T H n ) + t(d t H n ) -H n + n(n + a))(d t H n )(d t H n + N s ). (2.5) 

Remark 3. IfH„(T, t) is a function ofT only, i.e., H n (T, t) = Y n (T), then the PDE 
(2.3) reduces to the Jimbo-Miwa-Okamoto cr-form [29] associated with Painleve 
V: 

(TY n "f = (Y n - TY n ' + 2(Y n ') 2 + (v + v l +v 2 + v 3 )Y n 'f 
-4(v + *V)(vi + Y n ')(v 2 + Y n ')(v 3 + Y„'), 

where ' denotes differentiation with T, and 

v = 0, V] = -n- a, v 2 = -n, v 3 = -N s . 

Such a reduction can be obtained from a t — » oo limit in (1.20). 

Remark 4. A similar reduction can be found when H n (T, t) is a function oft only, 
i.e., H n (T, t) = Y n (t), which is a cr-form of Painleve V with parameters 

v = 0, vi=-n-a, v 2 = -n, v 3 = N s . 

This can be obtained from a T — » oo limit in ( 1.20). 



5 Dropping the T and t dependence notation for the sake of brevity. 
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With a change of variables the PDE, equation (2.3), can be converted into a 
form which will be convenient for the later computation of cumulants. Specifi- 
cally, let 

v 



T = - , (2.6) 

1 + cs 

t = v. (2.7) 

Under this transformation, H n (T, t) becomes H n (j^, v), which we write as H„(s, v) 
for the sake of brevity. Hence, 



H n (s,v) = vd v log M r (s,v). 
We have that v = t and s = -( j - 1) and consequently, 



d_ 

df 



(1 + cs) 2 d 



vc ds 
d_ (1 + cs) d d_ 

dt vc ds dv 



(2.8) 

(2.9) 
(2.10) 



Under the change of variables, the original PDE (2.3) becomes the following PDE 
in the variables (s, v): 



(1+csY 



vc 

(1+cs) 



(d s H n )(d v H n ) 



2(1 +csY 
(vc) 2 



(djlnf + (2n + N s + a + v){d v H n ) 



+ 2N s y -^^(d s H n ) - (2n - N s + a) K +CS \ d s H n ) ± A\{H n ) + A* 2 (H n ) 

vc v ' V v ' 

A\(H n )A* 2 (H n )(vc) 2 - (1 + cs) 2 [(vdl - d s )H n ][(l + cs)(vd 2 vs - d s )H„ + v 2 c(d 2 , v H n )\ 



2(vc) 2 [v(d v H n ) - H n + n(n + a)} 

where A\(H n ) and A* 2 (H n ) obtained from a corresponding transformation are 
(vc) 2 A\(H n ) 2 = (1 + cs) 4 [(vdi - d s )H n f 

+4(1 + cs) 2 ^v(d v H n ) -H n + n(n + a) 



(2.11) 



d s H n 



(I + csY(d s H n ) + N s vc 
(2.12) 



and 



(vc) 2 A* 2 (H n ) 2 = [(1 + cs)(vd 2 s - d s )H n + v 2 c(d 2 n ,H n )f 



+4^v(d v H n ) - H n + n(n + a) 
x ((1 + cs)d s + vcd^H n + N. 



((1 +cs)d s + vcd^H n 



vc 



(2.13) 
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We note that in the PDE (2.11), the only higher order derivatives are the mixed 
derivatives with respect to s and v, and second order derivatives with respect to v. 

It is of interest to note that the Hankel determinant in (1.20) or (1.23) has the 
following asymptotic forms: 

lim(r !A, 'D„[vv AF (-,r,r)]) = D n [w dLag (-,t,a,N s )], (2.14) 

\im(t-' ,N *D n [w AF (;T,t)]) = D n [w dLag (;T,a,-N s )], (2.15) 

(— »oo \ / 

where the generating weight is the "first-time" deformation of the Laguerre weight 

WdLagO, T,a,N s ) = x a e~\x + T) Ns . (2.16) 

In previous work by the authors [13], the Hankel determinant D n [wdhag(-, T, a, N S J\ 
was shown to arise in the analysis of the moment generating function of the single- 
user MIMO channel capacity, and was shown to involve a Painleve V differential 
equation. 

3. Coulomb Fluid Method For Large n Analysis 

In this section, we make use of the Coulomb Fluid method, which is particu- 
larly convenient when the size of the matrix n is large, to describe the MIMO-AF 
problem. The idea is to treat the eigenvalues as identically charged particles, with 
logarithmic repulsion, and held together by an external potential. When n, in this 
context the number of particles, is large, this assembly is regarded as a continuous 
fluid, first put forward by Dyson [22-24], where the eigenvalues were supported 
on the unit circle. For a detailed description of cases where the charged particles 
are supported on the line, see [21, 25, 42]. 

An extension of the methodology to the study of linear statistics, namely, the 
sum of functions of the eigenvalues of the form 

n 

Em)- 

7=1 

can be found in [21] and will be used extensively in this paper. The main benefit 
of this approach, based on singular integral equations, is that it leads to relatively 
simple expressions for characterizing our moment generating function. 

In Section 3.1, for completeness, we give a brief overview of the key elements 
of the Coulomb Fluid method, following [13, 21, 25, 42]. In Section 3.2, we 
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compute explicit solutions for the key quantities of interest in the Coulomb Fluid 
framework. Subsequently, in Section 3.3, we combine our Coulomb Fluid results 
with either ( 1 . 1 1 ) or ( 1 . 1 2) to directly yield analytical approximations for the SER 
of the MIMO-AF scheme under consideration. Quite remarkably, these are shown 
to be extremely accurate, even for very small dimensions. Similar accuracy is 
demonstrated in Section 3.4 for SNR cumulant approximations obtained from the 
Coulomb Fluid results. 



3.1. Preliminaries of the Coulomb Fluid Method 

We start by considering the ratio in the moment generating function expression 
(1.36), noting that it is of the form 



Z n(T\ f) _ e -[F H (T'S)-F n (fS)] 
Zn(f,t') 



where 



(3.1) 



1 C " " 

Z n (T', f) := - I exp -0(x u ...,*„)-£ /(*,-, T', t') ] ] dxj (3.2) 

[0,oo)' ! J 



n 

®(x u ...,x n ) := -2 ^ log \xj - x k \ + n ^ Vote) (3-3) 

\<j<k<n 1=1 

and F n := - logZ„ is known as the free energy. 

This expression embraces the moment generating function expression (1.36) 
with appropriate selection of the functions v (x) and f(x, T', f). A key motivation 
for writing our problem in this form is that it admits a very intuitive interpreta- 
tion in terms of statistical physics, as originally observed by Dyson (see [22-24]). 
Specifically, if the eigenvalues x\, . . . , x n are interpreted as the positions of n iden- 
tically charged particles, then 0(xi, . . . , x n ) is recognized as the total energy of the 
repelling charged particles, which are confined by the external potential nv Q (x). 
The function fix, T', f) acts as a perturbation to the system, resulting in a modi- 
fication to the external potential. The quantity F n (T', f) may be interpreted as the 
free energy of the system under an external perturbation f(x, T', t'), with F n (f, f) 
the free energy of the unperturbed system. 

For sufficiently large n, the system of particles, following Dyson, may be ap- 
proximated as a continuous fluid where techniques of macroscopic physics and 
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electrostatics can be applied. For large n we expect the external potential nv (x) 
to be strong enough to overcome the logarithmic repulsion between the particles 
(or eigenvalues), and hence the particles or fluid will be confined within a finite 
interval to be determined through a minimization process. For this continuous 
fluid, we introduce a macroscopic density cr(x)dx, referred to as the equilibrium 
density. Since v (x) is convex for xel, this density is supported on a single inter- 
val denoted by (a, b), to be determined later (see [42] for a detailed explanation). 
The equilibrium density is obtained by minimizing the free-energy functional: 

b b b 

F n (T',t'):= J o-(x)(n 2 v (x) + nf(x,T',t'))dx-n 2 J J cr(x) log \x - y\cr{y)dxdy, 

a a a 

(3.4) 

subject to 

b 

K)dx = 1. (3.5) 



j oix) 



With Frostman's Lemma [43, page 65], the minimizing cr(x)dx can be character- 
ized through the integral equation 

b 

n 2 \J {x) + nf{x,T',t')-2n 2 J log \x - y\cr(y)dy = A, (3.6) 

a 

where x e [a,b] and A is the Lagrange multiplier for the normalization condition 
(3.5), which can be interpreted as the chemical potential of the fluid [25, 42]. Not- 
ing that the integral equation above has a logarithmic kernel, taking a derivative 
with respect to x 6 (a, b) converts it into a singular integral equation of the form 

f'(x,T',n rcr(v) 

v (x) + = 2V dy, (3.7) 

n J x-y 

a 

where V denotes Cauchy principal value. 

Noting the structure (in n) of the left-hand side of (3.7), it is clear that cr(-) 
must take the general form: 

<r c (x,T',f) 

oix) = cr (x) + , (3.8) 
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where <To(x)dx is the density of the original system in the absence of any perturba- 
tion, while cr c (x, T', f) represents the deformation of cr (x) caused by f(x, T', t'). 
Furthermore, to satisfy (3.5), we have 

b b 
^ cr (x)dx = I, ^ cr c (x, T', t')dx = 0. (3.9) 

a a 

Substituting (3.8) into (3.7), and comparing orders of n, we see that <x (x) solves 

b 

Vo'Gc) = 2V f — dy, (3.10) 
J x-y 

a 



and cr c (x, T , t') solves 



b 

29 f P ' c(y,r,f/) rfy. 
J x-y 



Following [42], where the choice for the solution for cr has been extensively 
discussed based on the theory described in [44]; the solution to (3.10) subject to 
the boundary condition cr (a) = <x (£) = reads 



b 



V(fr -*)(*- a) f v' (x)-v' (y) 

o-q(x) = — = -dy, (3.12) 

2n J (x-y)J(b-y)(y-a) 

a 

together with a supplementary condition, 







y/(b - x)(x - a) 

a 

b 

The solution to (3.1 1) subject to j cr c (x, T', t')dx = is 



V (x) 

dx = 0. (3.13) 



b 



C J(b - y)(y - a) 

cr c (x,r,f) = - - . = P ^ '-f(y,T',1?)dy. (3.14) 

2n 2 V(5 - x)(x - a) J y-x 

a 
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Finally, the normalization condition (3.5) becomes 

b 



xv'Jx) 

-_dx = 2n. (3.15) 



V(£> - x)(x - a) 



The end points of the support of the density cr (x), a and b, are determined by 
(3.13) and (3.15), and will depend on parameters associated with v . For a de- 
scription see [42]. (We mention here a related problem, namely, the probability 
that there is a gap in the spectrum of the random matrix, has been studied using 
the Coulomb Fluid approach in [45, 46].) 

With the above results, for sufficiently large n, we may approximate the ratio 
(3.1) as (see [21] for more details) 

Z n (T',t') 



Z n {t',f) 



exp( - SfiT', - nSf(T', t')), (3.16) 



where 



b 



Sf(T',t') = I cr (x)f(x,T',t')dx, (3.17) 



b 

Sf(T',t') = -Jo- c (x,T',t')f(x,T',t')dx. (3.18) 

In the sequel, we will find explicit solutions for these quantities. 

3.2. Coulomb Fluid Calculations for the SNR Moment Generating Function 

Using (3.16), the moment generating function Aiy(T',t') in (1.36) takes the 
form 



M r (T',t') = - 



T'\ Mn Z n (T',f) 



r ) z n (t',o 



&xp\-Sf(T',t')-n 



Sf(T',f)-N s log\- 



, (3.19) 



for which, comparing (1.37) and (1.38) with (3.2) and (3.3), the functions v (x) 
and f(x, T', t') are identified as 

v (x) := x-fi\ogx, (3.20) 
f(x,T',t') := AUog&^j, (3.21) 
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where v (x) given by (3.20) is convex. 

From here, the Z-th cumulant can be extracted from the formula 



Ki = c 



(3.22) 

T'=t> 

The objective of the subsequent analysis is to evaluate the quantities Sf F and 
5^ F . As we shall see, we will need to solve numerous integrals which are not 
readily available. Thus, to aid the reader, we have compiled these integrals in 
Appendix B. 

We start by considering the end points of the support a and b. These are 
determined by equations (3.13) and (3.15), and satisfy 

ab = f, (3.23) 
a + b = 2(2+/3). (3.24) 

With the integral identities (B.1)-(B.3) in Appendix B, we obtain, after a few easy 
steps, 

a =2 + (3 - 2 Jl + fi, 

(3.25) 



b=2+/3 + 2 yi +j8. 

The limiting density ctq(x) in (3.12) can be computed using the integral identity 
(B.2) as 

V(fr - x)(x - a) 

cr (x) = , xe(a,b), (3.26) 

2nx 

which is the Marcenko-Pastur law [13, 47, 48]. Meanwhile, with the aid of the 
integral identity (B.4) in Appendix B, cr c (x, T' , t) given by (3.14) reads 

, ^ Ns ( W + a)(T'+b) W + a)(f + b) \ 

cr c (x, T ,t) = . 

2nyl(b-x)(x-a)\ x + T' x + f J 

(3.27) 

Using cr (x) and invoking the integral identities (B.5)-(B.7) in Appendix B, 
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gives 



Sf(7V) = 




N s ^[ab I (V^+ W+W+ty -t' 2 



t' -T + V(^' + a)(r' + ft) - V(f + a)(f + ft) 



( V^ft + W + «)(^' + ft)) - T' 2 



N s (a + b). I ( VC + a)(f + 55 + _ aft | N s (a + b), [V 




(y/(T' + a)(T' + b) + T') - ab) 4 \* ' 

(3.28) 



Moreover, using cr c (x, T , t) and invoking the integral identities (B.2), (B.8) and 
(B.12) in Appendix B gives 



3.3. SER Performance Measure Analysis Based on Coulomb Fluid 

Combining (3.19) with (3.25), (3.28), and (3.29) yields a closed-form asymp- 
totic expression for the moment generating function of the instantaneous SNR in 
(1.10). This, in turn, combined with either (1.1 1) or (1.12), directly yields analyt- 
ical approximations for the SER of the MIMO-AF scheme under consideration. 
The accuracy of these approximations is confirmed in Fig. 1, where they are com- 
pared with simulation results generated by numerically computing the exact SER 
via Monte Carlo methods. Different antenna configurations are shown, as repre- 
sented by the form (N s , Nr, Nd). The curves labeled "Coulomb (Exact SER)" were 
generated by substituting our Coulomb Fluid approximation into the exact expres- 
sion of the SER (1.11), whilst the curves labeled "Coulomb (Approx SER)" are 
generated by substituting our Coulomb Fluid approximation into the approximate 
expression for the SER (1.12). 

From these curves, it is remarkable that the Coulomb Fluid based approxima- 
tions, derived under the assumption of large n, very accurately predict the SER 
even for small values of n. This is evident, for example, by examining the set of 
curves corresponding to the configuration (2, 3, 2), for which n = 2. In fact, even 



4 V(r + a)(T' + b) W + aW + b) 



Sf(T',f) = -j-log 



( V<T' + a)(T' + b) + V0' + a)(f + b)f - (T - t) 2 j ' 

(3.29) 
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for the extreme cases with n = 1 , the Coulomb Fluid curves still yield quite high 
accuracy. 

We note that the results in Fig. 1 are representative of those presented previ- 
ously in [7, Fig. 5]. The key difference therein was that the analytic curves were 
based on substituting into either (1.11) or (1.12) an equivalent determinant form 
of the moment generating function to that given in (1.24) and (1.25). That rep- 
resentation, involving a determinant of a matrix with elements comprising sums 
of Kummer functions, is far more complicated than the Coulomb Fluid represen- 
tation, which is a simple algebraic equation involving only elementary functions. 
The fact that the Coulomb Fluid representation yields accurate approximations for 
the SER for small as well as large numbers of antennas makes it a useful analytical 
tool for studying the performance of arbitrary MIMO-AF systems. 




SNR (dB) 

Figure 1: Illustration of the SER versus average received SNR (at relay) y; comparison of anal- 
ysis and simulations. Each set of curves represents a specific antenna configuration of the form 
(N s ,N r ,Nd). For configurations with N s = 2, the full-rate Alamouti OSTBC is used (i.e., R = 1), 
whilst for N s — 4, a rate 1/2 orthogonal design is employed (i.e., R = 1/2). In all cases, QPSK 
digital modulation is assumed, such that M — 4. The relay power b is assumed to scale with y by 
setting b — y. 
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3.4. Coulomb Fluid Analysis of Large n Cumulants ofSNR 

In this subsection we use the Coulomb Fluid results to study the asymptotic 
cumulants of the SNR. We suppose J3 is fixed, and distinguish between two cases, 
p = and fi > 0. 

3.4.1. Case 1:13 = 

Here N R = N D . Then from (3.25), a = and b = 4, leading to 

Sf(T', t') - N s log |y J = ^{t'-T'+ Vr'(7" + 4) - VfC+4)) 



+2N s \og 



and 



Sf(TV') = ^-log 



7V + vr'C +4) J ' 
2 Vrr V(^' +4)(r + 4) 



(3.30) 



Vf 7 ? V(^' + 4)(f + 4) + 2r + 2r + rr 



.(3.31) 



From the closed-form approximation of the moment generating function (3.19) 
the cumulants will be extracted using (3.22). 
Recall that the variable T is related to s via 



r 



1 + cs 



where c = y/(RN s ). Therefore we expand log M y (f /(l +cs), f) obtained from the 
Coulomb Fluid formalism (3.30) and (3.31) in a Taylor series about s = 0. From 
(3.19), 

oo j 

logMJf/il +cs),f) = £(-1)^(0, (3-32) 



where the first few cumulants are 

cN : 



c 2 N s 

* 3 (Q 

c 3 N v 



- = 5(^ + 2- V^4f), 

——?—r + n(t' + l) - n Vr' 2 + 4r 1 1 —I, 

UN " = + n(3*' + 2) - n V?' 2 + 4?' 1 3 - - 
(f + 4) 3 \ f 



+ 4 (f + 4) 2 



(3.33) 
(3.34) 
(3.35) 
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K4(.Q 

c 4 N s 
c 5 N x 



and 

c 6 /V v 



174/V i / 3 

' + 4) 4 v ' \ f+4 



(f 

3120/V, 



+ 4 + 4) 2 (f + 4) 3 / ' 
(3.36) 



+ 12n(5*' + 2) 



(f + 4) 5 
-Yin yjt' 2 + 4f \5 - - 



8 



10 



+ 4 (t' + 4) 2 (t' + 4) 3 (t'+4) A 



(3.37) 



67440/Y, 



(f + 4) 



5 + 120n(3r' + 1) 



-120n V?' 2 + 4r\3 — 

\ f +4 



10 



14 



+ 4 (f + 4) 2 (f + 4) 3 (/"+4) 4 (f + 4) 5 /' 

(3.38) 

5.4.2. Case 2: fi >0 

In this case, the end points of the support, a and b, are given by (3.25). Hence, 
going through the same process, the first five cumulants are 



c 2 W, 



«f(f) 



- ^ + 2+0- ^(f +j8) 2 + 4f), 



((f + /3) 2 + 4f ) 2 



+ « 



(t' + \+^-n^t' +pf + 4t' 



(3.39) 

' 1 j3 z + (J3 + 2)f ' 
2 ((?' + jS) 2 + 4?') J ' 
(3.40) 



6N s (l+p)(f3 2 + (J3 + 2)t')t' 2 



+ n (3?' + p + 2) 

a 



((f + p) 2 + At'f 

I 1 / 2(1 +B+0 2 + (0 + 2)0 

-n J(f + p) 2 + 4f 3 - - , — 

2(0 i +p 2 + 2(/3+ l)(/3 + 2)?y 



c 4 N, 



((/' + /3) 2 + 4t'f 

36N/ 2 (\+p)((p 2 + fp+f)t' 2 + 2p 2 (p + 2)t'+p 4 ) 

((?' + j6) 2 + 4?') 4 

o I ^ — T/, 3f(f + 2+fi) 2t' 2 (t' -3-3/3) 

-3n J(t' + p) 2 + 4?' 1 + — — + ^ 

^ \ {f+p) 2 + 4t> {{f +p) 2 +4t'f 



(3.41) 
+ 3 n (4t' + p + 2) 
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2t' 3 (t' 2 + (3+p)t' -2-3p)\ 

(3.42) 



{{f + /3) 2 + At'f 



and 

K f(f) 24(W/ 2 (1 + fi)( (p 2 + ^p+^) t' 2 + 2/3 2 (2 + 0) f + p 4 ) ((2 + p)f + p 2 ) 

c 5 N s (( t > + pf + At') 5 

+ \2n(5t' +P + 2) 

-i2„J( f - + ^ + 4,-(i + 4 '' ( '"'; 4 ^- 8) 

> \ (t'+p) 2 + Af 

At' 2 (5t' 2 + 2(5p + A9)f - 37(0 + 1)) 

((f +P) 2 +At') 2 
2t' 3 [\0t' 2 + \2{fi + A)t' 2 - 3(125/3 + 334)?' + 348/3 + 232) 

((*' + py- + Aff 

10f' 4 (28f' 3 + 3(12/3 + 23)f' 2 - 2(135/3 + 242)?' + 158/3 + 79) 

+ -. 

((?' + P) 2 + At'f 

280m 5 (2t' 3 + (3/3 + A)t'{t' -A) + 5p + 2)\ 
((?' + P) 2 + At'f I 

The accuracy of these approximations is demonstrated in Fig. 2, where they are 
compared with simulation results generated by numerically computing the exact 
cumulants via Monte Carlo methods. As before, different antenna configurations 
are shown, as represented by the form (N S ,N R ,N D ). From these curves, the ac- 
curacy of our Coulomb Fluid based approximations is quite remarkable, even for 
small values of n. 

We note that exact expressions for K\ and k 2 (i.e., the mean and variance) were 
derived previously in [7, Corollaries 1 and 2], and these were expressed in terms 
of summations of determinants (for the mean) as well as a rank-3 tensor (for the 
variance), each involving Kummer functions. Such results are obviously far more 
complicated than the Coulomb fluid based cumulant expressions in (3.39)-(3.41), 
which involve just very simple algebraic functions. 
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Figure 2: Illustration of the mean K\, variance Ki, and third cumulant /C3 of the received SNR at the 
destination y, as a function of the average received SNR at the relay y. In each case, the simulated 
cumulants are compared with those obtained based on the Coulomb fluid approximation. As in 
Fig. 1, each set of curves represents a specific antenna configuration of the form (N s , Nr, Nq)- For 
configurations with N s = 2, the full-rate Alamouti OSTBC is used (i.e., R = 1), whilst for N s = 4, 
a rate 1 /2 orthogonal design is employed (i.e., R = 1/2). The relay power b is assumed to scale 
with y by setting b = f. 



4. Power Series Expansion For H n (s, v) 

In this section we seek power series solutions for (2.11), to determine the 
cumulants. We begin by removing the square roots of the PDE and arrive at the 
following lemma. 

Lemma 1. Let 

K := 2(v(d v H n )-H n + n(n + a)), (4.1) 

and 

(l+cs) 2 , tt w tt v 2(1 +cs) 3 



L := K 



H n - 2'-^L(d s H n )(3 v H n ) - ^^(d s H n y 
vc v /v ' (vc) 1 v ' 



(1 + cs) s(\ + cs) 
-(In + N s + a + v)d v H n - 2N S - -d s H n + (In - N s + a)— -d s H n 

VC V 



26 



~ (1 (vc" )2 [ (v ^ J ~ ds)H "}[ {l + cs) ( v5 v.v - W + v 2 c{d 2 vv H n )\ (4.2) 
then the PDE (2,11) may be rewritten in an equivalent square-root-free form, as 

AK 2 A\(L + A\f = (L 2 + K 2 (A\ - A\) - A 2 Aj) , (4.3) 
where A i andA 2 are given by (2.12) and (2.13) respectively. 

Now, recall from (2.8) that the function H n (s, v) is related to the moment gen- 
erating function (1.20) through 

H n (s,v) = vd v \og M 7 (s,v), (4.4) 

and note that the logarithm of the moment generating function has the expansion, 

7=1 J ' 

where Kj(v) is the jth cumulant. Note that for the sake of brevity we write M r (s, v) 
in place of M y [r^, vj. Consequently, from (4.4), H n (s, v) has the following ex- 
pansion in s, 

fl' f (v) 

H n (s, v) = v 



^(_|| V/ \/^ r .v'. (4.6) 

7=1 J ' 



where ' denotes differentiation with respect to v, and a 7 (v) is related to the 7th 
cumulant Kj(v) by 6 

a/v)c y iV, = Kj (v). (4.7) 

4.7. Analysis of K\ 

In this subsection, the computation of the mean K\ is presented. To this end, 
upon substituting (4.6) into the PDE (4.3), keeping the lowest power of s in the 
resulting expansion gives a highly nonlinear ODE expressed in the factored form 

Tj 1 ^ = 0, (4.8) 



6 We introduce aj in order to make a clearer comparison to the Coulomb Fluid model and [7] 
without having to exactly specify a value for the constant c, since it falls out of the subsequent 
expansion. 
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where 



:= (a'Hv)) 2 (a'Hv)) 2 + n(n + cx)In s 2 (a\(v)f - N s 2 (a\(v)) - n(n + a)) 



(4.9) 



and 



T 2 



; (<'(v)) + v(oi'(v)) + 4n(n + ^(v) - | 



-(2n + v + a)' 



v 2 (a;'(v)) 2 + 4n(n + a) (^(v)) (^(v) - 



(4.10) 



Hence, (4.8) is equivalent to either = or = 0. 



The ODE W] = has two solutions, the first of the form a\(v) = c\v + C2 
where c\ and C2 are constants to be determined. The second is more complicated 
but can be succinctly written as 



ai(v) = f{c\ , N s , n, a) v sin (2 V« V" + ce log vj 

+g{c\,N s ,n,a) v cos (l^Jn^ln + a log v) + C2 + (4.11) 



where / and g are the integration constants. 

Whilst these are valid mathematical solutions, it is clear that by taking n large, 
they differ drastically from that predicted by the Coulomb Fluid method and, for 
small n, they differ with the solutions obtained in [7]. This suggests that these are 
not the solutions of interest to our problem; thus, we set them aside and henceforth 
focus on ¥^ = 0. 

It becomes evident at this point that we require initial conditions for the cumu- 
lants Kj(v)/(c j N s ), since the PDE (4.3) gives rise to a system of ODEs satisfied by 
Kj(v)/(c^N s ). However, (4.3) is not supplemented by any initial conditions from 
which one deduces the initial conditions at f = of the system of ODEs. 

As we shall see, results from the Coulomb Fluid analysis (3.39)-(3.41) are cru- 
cial, as these will specify initial conditions at f = 0. Consequently, the Coulomb 
Fluid results are in fact leading order approximations to the exact results. With- 
out the Coulomb Fluid analysis, each cumulant Kj(v)/(c j N s ) would carry unknown 
constants of integration. 

We now examine the equation Wi' = in the large n region, which leads to an 
asymptotic characterization of a\{y) and thus a^. To proceed, we scale the variable 
v by v = nt' . Also, note that a = n(m/n - 1) = nfi. We find that ai(f) satisfies the 
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following ODE: 
= 



t' 2 t' 12 

— (a\"(f)) + - «(?')) + 2/i(l + j8) (2a; (f') - /i) 
w n 



-(f + 2+/5Y 



t' 2 (a'lit')) 2 + An\\ +{3) (a[(f)) (a[(f) - n) 



For large n, keeping the leading order term, we arrive at 

i2 ,2/,/,o , o\2 



da\ n 
1F~2 



n\t' + 2+ /3Y 
A{f +p) 2 + \6f' 



(4.12) 



(4.13) 



whose solutions are 



ai (f) = Ut'+ V(?'+y8) 2 + 4f) + C 1 



where C\ is a constant of integration. We retain the solution that is bounded as 

t' -> oo, 



Cl\ 



(4.14) 



Comparing this with the corresponding result (3.39) obtained from the Coulomb 
fluid approximation, we see that k^ f (0) = cN s n, which suggests ai(0) = n and 
hence C\ = n(l + fi/2). Therefore, the first cumulant K\ to 0(n) becomes 



Large n (f) 



cN s 



= 

= n -{2+fl + t'- V(f +j8) 2 +4r). 



(4.15) 



4.2. Analysis of Ki 

Having characterized K\, we now turn to the variance k 2 . To this end, equating 
the coefficients of the next lowest power of s to zero in the expansion of the PDE 
(4.3), results in an ODE which can be factored into the following form 



(4.16) 



The ODE = 0, makes a reappearance, which we set aside, leaving the ODE 



v P (2) + v P (2) +¥ (2) 



o, 



(4.17) 
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where Wf\ Wf } and are given by 

8v 2 \ct\ - ^{v 2 ( 2 « - v + a) (ai") 2 + (2n + v + a)n{n + a)a\ (a\ - l))fli' 



2A^ 



r 



(2) 



-v 4 (2« + v + or) (fli") 4 + 2v 3 (2n - v + or) [ai - - ) (ai") 3 



-v 



((2n + v + a)((a + v)(a - v) + 20n(n + a)) - 16n(n + a)(2n + a)) (oi') (a\ - 1) 



— n (n + a) (2« — 3v + a) 



K') 2 



+8 (2« + v + a)n(n + a)a\ (a\ — 1) 

ai' - ijfli" + ((v + a) 2 + 4v«)ai'(l - ai') + « (« + a) 



(4.18) 



2v 4 (2n(« + or) — 1) (a\") 2 - v 3 (va2 + 8n(n + a) (02' — l))ai" 
+4v 2 n(« + a)(v|«i' - ^j«2" + ((v + a)(2« + v + a) + 2«(4« + a - - 1) 

+n(n + a) + (a 2 ' - l)ai' - ^ a 2'j a\" + 4 v 3 ^«(« + ar) - i j (fli") 3 - v 4 a"' (a'/) 2 

v (2n + v + a) 2 02" - 2^4(2« + v + a) 1 - 8n(« + a))n(« + a) - v(2n + a) - a 2 jai' 

(«i") 2 



-(l2«(« + a) + 2vn + (v + a) a)a 2 ' + 4n (n + a) ((2n + v + a) 2 - 2n(« + a) + 3 J 
-4«(« + ar)vV/ la', - i J (vaj" - a 2 ' ) ~ 12vn 2 (« + a) 2 a", 



(4.19) 



and 

vp(2) 
4n(n + a) 



((2« + v + a-) 2 + 4n(n + a) - 2) {a x ') 2 - ((2« + v + a-) 2 - 4n(« + a) + ^)a 2 ' 
+M 2(2« + v + a) 2 + 1 - 8n(n + ar))a 2 ' - 3(2n + v + a) 2 + 1 + 4n(n + a)^ai 

+v [Anv + (v + a) 2 ) (ai') (1 - a[) + n(n + a) a 2 " 
+(2nv + (v + a) ff)a2'(l — a' x )a\ 

2(2n + v + a){2n + a — 2(2n + v + a)n(n + a)) 



a i 



+8«(« + a)(2n{n + a) - 1) 



(fli') 2 (a'i - 1) 
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+2((2n + v + a) 2 - An(n + a))n(n + a) («i') 2 - 2n (n + or) |- i« 2 ' + » (« + 
+l(6n(n + a) - (In + v + a-) 2 - 1 )n (n + a) a i ' . (4.20) 

Clearly, ^Fj depends on a\ only, whilst and W\ depend on both a\ and £?2. 
To extract the large n behavior of K2, we replace v by nt' in (4.17) and make use 
of the large n formula for a\(t') in (4.15), resulting to a highly non-linear ODE 
satisfied by a2(t'). To proceed further, keeping the highest powers of n, a first 
order equation is obtained for ci2(f), 



da%{f) 
df 



- n J(t' + f3) 2 + At' 



t' +13 + 2 



t'(fi+\) 



(t'+/3) 2 + Af ((f + f3) 2 + At') 2 



(4.21) 



((?' + (3) 2 + At') 3 

Integrating this leads to an expression for a2(t') which, again yields an unknown 
integration constant. This constant can be determined through a comparison with 
the Coulomb Fluid results for K 2 (t') in (3.40) at f = 0, giving a 2 (0) = n. 

The above analysis gives the leading order characterization of the variance, 



= n (t' + 1 + I ) - n yfft' + /3) 2 + At' 



i 1 /3 2 + (J3 + 2)t' 
2 ((?' + P) 2 + At') 



N s (l +{3)t' 2 



((*' + f3) 2 + At') 



,\ 2 ' 



(4.22) 



4.3. Beyond K\ and k 2 

The same procedure for computing k x and k 2 easily extends to higher cumu- 
lants. Here we give the leading order formula for /c 3 , /c 4 and k 5 , 

= a 3 (0 



-n J(t' + (3) 2 + At' 



x 



_ 2(\+/3+f + (J3 + 2)t') + 2{p 3 +J3 2 + 2(J3+ 1)(J3 + 2)0 



+n (3f' +/3 + 2) + 



(t'+J3) 2 + Af {{f + f3) 2 + At'f 

6N s t' 2 (l + /3)(J3 2 + (J3 + 2)t') 



((?' +p) 2 +At'y 



(4.23) 
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Large n 



it') 



c 4 N s 



and 



^ argen (f') 



c 5 N x 



an(t') 

-3n yj(t' + /3) 2 + 4 f / 



x 



W + 2+/?) + 2?' 2 (?' - 3 - 3/3) 2?' 3 (?' 2 + (3 + /3)?' -2-3/3) 



if + + 4?' ((f + £,2 + 4f / f ( (f + £,2 + Af f 

36N s t' 2 (1 + /3)((/3 2 + + 2 |)f' 2 + 2/3 2 (/3 + 2)?' + 



+3«(4?' + /3 + 2) + 



((?' + /3) 2 + 4?') 4 



(4.24) 



05(0 

-12n J(? + /3) 2 + 4? 



1 + 



4^(11^-48-8) 



+ 



(f + yS) 2 + 4?' 

At' 2 (5?' 2 + 2(5/5 + A9)t' - 37(J3 + 1)) 

((/' + /3) 2 + 4?') 2 
2t' 3 (lO?' 2 + 1208 + A)t' 2 - 3(125/3 + 334)?' + 348/3 + 232) 

((*' + /3) 2 + 4?') 3 
10?' 4 (28?' 3 + 3(12/3 + 23)?' 2 - 2(135/3 + 242)?' + 158/3 + 79) 

((?' + /3) 2 + 4?') 4 
280w 5 (2?' 3 + (3/3 + 4)?'(?' - 4) + 5/3 + 2) 



+ 12«(5?' +/3 + 2) 



((?' + /3) 2 + 4?') 5 

24(W 4 ?' 2 ( 1 + fi){ (0 1 + ^3 + t' 2 + 2/3 2 (2 + fi) ?' + /S 4 ) ((2 + /S) ?' + /3 2 ) 



((?' + /3) 2 + 4?') 



/\5 



(4.25) 



4.4. Comparison of Cumulants obtained from ODEs with those Obtained from 
Determinant Representation 

This subsection serves as a check for consistency of our equations. For small 
values of n we compute the Hankel determinant from the moments formula ( 1 .25), 
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since D n {T,t) = det(ji i+ j(T,t))" jl . The moment generating function in s and t 
reads 



det(c i+ ,O,0) 



n-l 



My(S, t) = 



( 



1 + CS 



1 




det(c I+; (0, 0) 




(4.26) 



where 




(4.27) 



which was derived in [7]. For small fixed integer values of n and N s , e.g., n = 2, 3 
and N s = 10, and n = 4, 5 and Af v = 1, the above determinant can computed 
without much difficulty, from which the cumulants follow. It can be seen that 
^ = ai(t) and ^ = a 2 (t) obtained from (4.26), satisfied third order ODEs for 
a\{t) given by (4.10) and a 2 (t) given by (4.17). Similar results hold for the higher 
cumulants, which provide a consistency check. 

5. Large n Corrections of Cumulants obtained from Coulomb Fluid 

We have shown that the PDE (4.3) satisfied by 



may be used to generate a series of non-linear ODEs that are satisfied by the 
cumulants K/(v). Under a large n assumption, where v = nf, the first few of 
these ODEs are approximated as first order ODEs for A^ arge "(?'), whose solutions 
matched exactly with that obtained from the Coulomb Fluid analysis for the cu- 
mulants Kf F (t'). 

We give a method in this section where the non-linear ODEs generated from 
the PDE (4.3) are employed systematically to obtain "correction terms" to the 
Coulomb Fluid results. 



H n (s,v) = vd v \og M y (s,v) 
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5.1. Large n expansion of K\ 

We assume the first cumulant or K\ has the following large n expansion 



ai(t') = 



cN s 

' e k {f) 



= ne^(f) + e (f) + y^Y-. (5.1) 

U n 

Substituting the above into (4.12), and setting the coefficients of n j to zero, a 
system of first order ODEs for ek(t') is obtained. These are are solved successively 
starting from e-\(t'), followed by en{t') and so on. 

The expression for /^ aige n (f) given by (3.39), gives rise to the initial conditions 

e_i(0) = l and e k (0) = for * = 0,1,2,3,.... (5.2) 

As a result e-\(t') is found to satisfy the ODE 

<de. l {f) n\ 2 n 2 (t' + 2+/3) 2 



df 2) 4(t'+p) 2 + l6f 



and has two solutions 



e-tf) = l -(2+(3 + t'± ^/(t'+/3) 2 + 4t') , 



both satisfying the initial condition ei(0) = 1. 
We retain the solution that is finite at infinity 



(5.4) 



to match with that obtained from the Coulomb Fluid computation, namely, (3.39). 
Continuing, we set the coefficient of n 3 to zero, which implies, 

den 

-£ - °- <«) 

Setting the coefficient of n 2 to zero gives rise to an ODE involving e_ 1 (?')> e '-\{t'), 
e"l\(f \ e ^') anc ^ Simplifying, with given by (5.4) and e' Q (t') given 

by (5.5), we find that ei(t') satisfies: 

((t' +B) 2 + 4t') 7 ' 2 de x , 2 , 
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Setting the coefficient of the next lowest power of n to zero gives rise to an ODE 
involving e_ x (f), e'_ x (t'), e'"(f), e' (t'), e' Q (t'), e" (t'), e\(f), and e' 2 {t'). From the 
expression of and e' (f) given by (5.4) and (5.5) respectively, we find 

dei 

Continue with this process, the next three terms e 3 (t'), e 4 (t') and e 5 (t') are found 
to satisfy the following equations: 



{(? +p) 2 +4t') i3/2 de 3 
t'il+P) df 

IF 

((t'+/3) 2 +4t') 19/2 de 5 
t'{l+P) df 



l?\t'), (5.8) 



0, (5.9) 
4 5 V), (5.10) 



where l[ (?) and l\\t') are given by 

f$(f) = -40t' 6 - 16(6 + 2)t' 5 + (149/3 2 - 79/3 - 79)?' 4 + (fi + 2)(145/3 2 - 27/3 - 27)?' 3 
-(19/3 2 - 236/3 - 236)/3 2 t' 2 + 37(/3 + 2)/3t' - 2/3 6 , (5. 1 1) 

if V) = - 1260?' 10 + 412 (/3 + 2) f' 9 + 12 (325/3 2 - 97/3 - 97) f' 8 

+806 + 2)(1999^ 2 - 516y3 - 516)f' 7 

- (7967 Z3 4 - 70762/3 3 - 61492/3 2 + 18540/3 + 9270) t' 6 

-(/3 + 2)(23233/3 4 - 4375Q8 3 - 41500/3 2 + 4500/3 + 2250)f' 5 
-2(4294^ + 33485 /3 3 + 10575^ 2 - 45820/3 - 22910)jSV 4 

+34(/3 + 2)(107/3 2 - 695/3 - 695)^ t' 3 + (1717/3 2 + 10072/3 + 10072)/3 6 ?' 2 

-217(/3 + 2)/3V + 2/3 10 . (5.12) 

Solving ODEs (5.5)-(5.10) with initial conditions e,t(0) = 0, we obtain, 

Ks(t') K^it') 9 / : — /i\ 

S£ . ^2 +{1+ ^^ + ^ +4f ^^2 + o{-), (5.13) 
where Afty), A^ 1 V) and A<V) are given by 

A<V) = ; 7, (5-14) 

((?' + /3) 2 + 4f ) 
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A (V) = 1 , | W 2 (l +/ ?) 

3 ((f + y6) 2 + 4f') 4 ((*' + /3) 2 + 4f') 5 ((*' + pj 2 + 4f') 6 ' 



4V) = 7 + 



2t'(331t' - 3Sp - 76) 



((?' + y6) 2 + 4f ) 5 ((f + p) 2 + At') 6 
33f' 2 (15f' 2 + 484f' + 60t'p - 106/3 - 106) 
((*' + p) 2 + At') 1 

2002?' 3 (6?' 2 + 52?' + \5t'p - 18/? - 12) 50050f' 4 (f' 2 + At' + 2t'p -2/3-1) 
((*' + p) 2 + At'f {{t' + p) 2 + At') 9 

(5.16) 



5.2. Large n Expansion of 

In computing a series expansion for the variance, we proceed in a similar way 
as was just done for the mean. First, we substitute the expansion for a\(t') from 
(5.13) and 



a 2 (f) 



*2(Q 

»/-i(0 + /o(0 + y=*Hr, 

h n 



(5.17) 



into (4.17). Equating the coefficients of n k to zero, a system of first order ODEs 
for fk{t') are found. Comparing with the Coulomb Fluid results for K2(f) in (3.40) 
gives rise to the initial conditions 

/_i(0) = l and /*(0) = for jfe = 0,1,2,3,.... (5.18) 

In this case, setting the coefficient of n 12 to results in a first order ODE for f-i(t'): 



df_ 

df 



1 - 



f + J3 + 2 



2Q3 + \)f 



X(f + P) 2 + At') 111 ((f + P) 2 + At'f 12 
The solution of this, with the initial condition /-i(0) = 1, reads 

( 



(5.19) 



f-i(f) = ( f/ + 1+ f)- ^(t'+I3) 2 +Af 



1 J3 1 + (fi + 2)f 

2 ((f + p) 2 + At') 



(5.20) 
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Continuing the process, we obtain an ODE for fait'): 

df 2N s (J3+l)t'(t'-/3)(t'+/3) 



df 



((f + p) 2 + At') 



3 



(5.21) 



The solution with the condition /o(0) = reads 



. ^ (522) 

((f + £) 2 + 4r) 2 



from which it can be immediately seen that 

*/-i(0 + /o(0 = -^-, 

and we recover i^(f) found previously. 

The ODEs satisfied by fk(t'), k = 1, 2, 3, 4 are reported in Appendix C, equa- 
tions (C.1)-(C.4). These will allow us to compute fj(f) with appropriate initial 
conditions. 

Briefly, the idea is that to determine the ODE satisfied by the k-th correc- 
tion term fk(f), for k odd, the preceding j(k + 1) ODEs satisfied by fj(f), j = 
-1, 1, 3 . . . , k are employed. For even k, the previous | ODEs satisfied by fj(t'), 
j = 0, 2, 4 ... ,k are employed. 

Going through the procedure described, we find that K 2 (t') has the following 
large n expansion, 

Kj (f) K^ F (t') , / 2 A (2) ( ^, ) 

+ (l+ y Sy 2 > /a'+yS) 2 + 4?'2] 



c 2 W, c 2 N s v " > v " £^ n 2 *" 1 

2 D (2) 



+ ^ 1+ ^' 2 Z^ + °b • (5-23) 

&=1 



where 



A?V) •= ? 5t'(t' + 2 +P ) 

1 ((f + ^) 2 + At') 3 ((f + ^) 2 + At') 4 ' 

:= ! W-2/3-4) 104^(1 + 

((*' +/3) 2 + 4f) 3 ((*' + /3) 2 + 4f) ((?' + /3) 2 + 4?') 
A ( 2 V) .= 3 , 7^- 9/3 -18) 

3 ' ((f +p) 2 + At') 4 ((t' + ft) 2 + At') 5 
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2lt' 2 (9t' 2 + 13? + 9t'/3 - 49/3 - 49) 1155?' V + 3? + t'p -3/3-2) 
((f + p) 2 + At') 6 ((f + p) 2 + At') 1 

(5.26) 

2*' (337*' - 38/3 - 76) 



Bf\f) := ^—7 + 



((?' + p) 2 + At') 4 ((?' + p) 2 + At') 5 
t' 2 (A95t' 2 + \59AAt' + 1980?'/3 - 3496/3 - 3496) 

((? + P) 2 + At') 6 
56t' 3 (2lAt' 2 + 1853?' + 535?'/3 - 642/3 - 428) 

((*' + p) 2 + At') 1 
49840?' V 2 + At' + 2t'p - 2/3 - 1) 

— tL r JL — • ( 5 - 2? ) 

((*' + p) 2 + At'f 



5.3. Beyond K\ and k 2 

The procedure adopted above for computing the large n expansion series for 
K\(t') and K 2 (t') easily extends to the higher cumulants. By way of example, here 
we focus on the third cumulant k 3 . In this case, an asymptotic expansion for a^(f) 
is assumed, 

a3i0 = ^ 

= ng^iO + goiO + Y^, (5.28) 

m n 

along with the initial conditions, 

g-i(0) = 2 and g k (0) = for Jk = 0,1,2,3, (5.29) 

In this case, we find that the large n expansion of the third cumulant reads 

„CFf t '\ , 2 ,(3) 



+ (l+P)t' 2 yjif +Pf+At'Y j ' 

2 o(3)/ 



c 3 N s c 3 N s 



2k-l 

k=l 



k=\ 
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where 



12 10f'(f' + 4/3 + 8) 140f' 2 (l+/8) 



1 ((f + P) 2 + At'f ((f + /3) 2 + 4f') 4 ((f + /3) 2 + 4f') 3 

J 2f'(f'-2/3-4) | 32^(1 + fl \ 
S I ((f + ^) 2 + 4f ') 4 ((f + y8) 2 + 4f') 5 / ' 
B (3) (f ') = 6 | 6f(37f -191-38) 

2 ((f + p) 2 + 4/') 3 ((f + /3) 2 + 4f') 4 

12f' 2 (21f' 2 + 172*' + 2 If'/? - 134/? - 134) 1560f' 3 (f' 2 + 3f' + f£ - 3/? - 2) 
((f + P) 2 + 4f') 5 ((f + P) 2 + 4t'f 

(5.32) 

12 42f'(35f' -8/3- 16) 126f' 2 (llf' 2 + 206f' + 26f'yS - 67/3 - 67) 



A (3) (f') = - 

3 ((f + /3) 2 + 4f') 4 ' ((f + p) 2 + 4f') 5 ((f +/3) 2 +4f') 6 

924f' 3 (21f' 2 + 152f' + 45f'/3 - 63/8 - 42) 60060f' 4 (f' 2 + 4f' + 2f'/3 - 2/3 - 1) 
((f + /3) 2 + 4f') 7 ((/' + P) 2 + 4t'f 

2 /20f'(12f' - p - 2) 36f' 2 (5f' 2 + 166f' + 20f'/3 - 34/3 - 34) 
I ((/' + /3) 2 + 4f') 5 ((f + j8) 2 + 4f') 6 

12f' 3 (378f' 2 + 3317f' +945f'/3- 1134/3- 756) 
((f + P) 2 + At') 1 



+N 



19392f' 4 (f' 2 + 4f' + 2f'yS -2/3-1) 



((f + y3) 2 + 4f) 



s 



(5.33) 



g4 3) (fQ 1 f'(1199f' -97^- 194) 

6 ((f + p) 2 + At'f ((f + P) 2 + 4f) 5 

f' 2 (3074f' 2 - 42140f - 5340f'/3 + 6004/3 + 6004) 
((f'+/3) 2 +4f) 6 

f ' 3 (4455f' 3 + 40366f' 2 - 460496f' + 4455f' 2 /3 - 130120f'/? + 94380/3 + 62920) 

((f'+/3) 2 +4f) 7 

28f' 4 (2247f' 3 + 4309f' 2 - 71182?' + 2247f' 2 /3 - 32965f'/3 + 19104/3 + 9552) 

((f'+/?) 2 +4f) 8 
_ 199360f' 5 (f' 3 - 15f' + t' 2 p - \0t'p + 5/3 + 2) 
((f + /3) 2 + 4f') 9 

For ease of reference, the ODEs satisfied by gk(t'), k = 1,2, 3,4, are placed in 
Appendix C, equations (C.9)-(C. 12) 

In summary, this entire section has shown that by expanding the cumulants 
K\(t'), K 2 (t') and K 3 (t') into asymptotic series in n, the Coulomb Fluid results are 
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recovered as the leading order contributions in the large n scenario. It is also seen 
that, by examining the finite-/? correction terms for each cumulant, no terms of 
0(n 2 ) or higher are present within the expansions. 



6. Asymptotic Performance Analysis Based on Coulomb Fluid 

In this section, we return to the analysis of the moment generating function, 
and consider the high SNR scenario (i.e., as y — » <x>). To this end, we will study 
the Coulomb Fluid based approximation derived in Section 3 (to be complemented 
in the following section through analysis based on Painleve equations), where the 
variables T and f are taken to be dependent on y, namely, 

T(s,y) = - t{y) = z , where b := vy . (6.1) 

1 + — n b 

Note that as y — » 00, 

N R 
nv 

To obtain the desired high SNR expansion, we compute the moment generating 
function M y as s — » 00 and y — > 00. The Coulomb Fluid based representation 
(3.19), when expressed in terms of T'{s, y) and — reads 



Myis) « e X p(- 5 f(rv,n^)-4^Hr(,,y),^)-^iog(" yr( ^ v 

\ v nv' v nv' \ 



Nr 



(6.2) 



where Sf F and 5^ are given by (3.28) and (3.29) respectively. Our goal is to 
compute the expansion of M y (s) as s — » 00 and 7 — > 00, which turns out to 
be an expansion in (ys) 1 . It turns out that the cases fi = and + behave 
fundamentally differently, and as such, these are treated separately. 

6.1. The Case of/5 = 

With /3 = 0, we have Nr = n. An easy computation shows that M y admits the 
following expansion: 

00 

M r^ = Z (^2- (6 - 3) 
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where A , A\, A 2 , A 3 , ... are constants independent of s and y. The leading expo- 
nent d is given by 



d = N, 



(6.4) 



whilst the first few Ar are 



= — J ^vexp 1- VTT4^ , (6.5) 



■A 2 = 



N s (l +4v) + 8«(2v- 1) 



RN< 



-A, = -4- 



2N S (n, Vl +4v - 4«) ■ 

v 7 [N s VTT4T - 4«) 

N. s (n s Vl +4v - 4«) - Vl +4v 

[32jV. 5 n 2 - Sn(2N s 2 + l)(2v - 1) - 3iV s (l + 4v)] VTT4T 
2N S (N s Vl +4v - 4n) 2 - iV s (l + 4v) - 8n(2v - 1) 
32«[^(1 + 4v) - 2«JV s (2v - 1) - 3v + ±] 



(6.6) 
(6.7) 



(6.8) 



2N S (N s VTT% - 4n) - N s (l + 4v) - 8«(2v - 1) 

We have refrained from presenting A k , k > 4, as these are rather long. 

Remark 5. For the special case b = y or v = 1, implying equal relay and source 
power, Aq reduces to the remarkably simple formula: 



(RN s ) N *( n -^(p 2N ' n lN s n 
A = — exp 



20^ 

where ip = (1 + V5)/2 is the Golden ratio. 



(6.9) 



6.1.1. High SNR Analysis of the Symbol Error Rate ( SER ) 

Based on (1 . 1 1), the SER of MPSK modulation can be expanded at high SNR 
using (6.3), resulting in 



MPSK 



I oo 



& (rsMPs K y' +f/2 



(6.10) 
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where 





i d A®) = j $™ ld+t ode . (6.11) 

o 

Considering the first order expansion, following [49] we may write 

/Wsk = (G a y) G " + o(y~ Gd ), (6.12) 

where we identify 

G d = N s (n-^\ (6.13) 
as the so-called diversity order, and identify the factor 

G a = gMPSK I I (6. 14) 

as the so-called array gain (or coding gain). We note that the result for G d above 
is consistent with a previous result obtained via a different method in [8], whilst 
the expression for G u appears new. 

Whilst it appears that a closed-form solution for the integral (6.11) is not forth- 
coming in general (though it can be easily evaluated numerically), such a solution 
does exist for the important special case of BPSK modulation, for which M = 2. 
In this case we have the particularization, = n/2, for which [50] gives 

W) = T nd + m + D ■ (6 - 15) 

Hence, using (6.14), G a admits the simplified form 

The high SNR results above are illustrated in Fig. 3. The "Simulation" curves 
are based on numerically evaluating the exact SER relation (1.11); the "Coulomb 
Fluid (Exact)" curves are based on substituting (6.2) into (1.11) and numerically 
evaluating the resulting integral; the "Coulomb Fluid (Leading term only)" curves 
are based on (6.12); whilst "Coulomb Fluid (Leading 4 terms)" curves are based 
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on the first four terms of (6.10). The leading-order approximation is shown to 
give a reasonably good approximation at high SNR, whilst the additional accuracy 
obtained by including a few correction terms is also clearly evident. 




SNR (dB) 



Figure 3: Illustration of the SER versus average received SNR (at relay) y; comparison of analysis 
and simulations. Results are shown for Nr = No = n and N s = 2, with the full-rate Alamouti 
OSTBC code (i.e., R = 1). QPSK digital modulation is assumed, such that M = 4. The relay 
power b is assumed to scale with y by setting b = |y. 



6.1.2. High SNR Analysis of the Probability Density Function ofy 

With the moment generating function expansion given above in (6.3), we may 
also readily obtain an approximation for the probability density function (PDF) of 
y, denoted f r (x), by direct Laplace Transform inversion. In particular, we obtain 



f=0 

For the case of very large y, with 



fyW = L T{ d + (i2) ■ <6 ' 17) 
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the leading term gives an approximation for the PDF deep in the left-hand tail. 
Of course, with the inclusion of more terms, a more refined approximation is 
obtained. 



6.2. The Case of/3 ± (with N R < N D ) 

For the situation where 0, two sub-cases arise. This first is N R < No, for 
which N R = n; the second is N R > N D , for which N R = (1 + f3)n. 

In the following, we will focus on the first sub-case, N R < N D . It turns out 
however, that our results also apply for N R > N D upon transforming the quantity 
v to v* by 



(6.19) 



1+/*' 



where v*(> 0) is interpreted as the (fixed) scaling factor between b and y, i.e., 

(6.20) 



b 

v* = - . 

7 



The moment generating function given by (6.2) admits an expansion distinct 
from the/? = case which does not have fractional powers of 1/Cys), reading 



AL(s) 



DO , 



ci+l 



(6.21) 



/=o 



where 



d 



nN< 



(6.22) 



and 



* (jv —nB) I — 

(1 +v/?) 2 +4v) 2 " (l +2y + vy6+ y(l + y/3) 2 +4y) 2 



^(2+/3) 



2 ^(2n+N s ) 



Ai = 



2v/3 2 



y'^((l + y/3) 2 + 4y)~(l + /?)^ j8 
exp(-^(l + vjS - 7(l + v/3) 2 + 4y)) , 

iy S ^(1 + y^) 2 + 4y - (2n + AQ0(1 + vfi) - 2N S 



(6.23) 
(6.24) 
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In this situation, the sub-leading terms are very complicated, however, the jth 
term in the expansion can be written in the following form: 



A/ = 



A Q R j N s j+i 
2(j!)(v^ 2 ) 



- [EjN s fi Vd+v>S) 2 + 4v + Fj] , (6.25) 



where Ej and Fj also depend upon v,N s ,n and /?. 
In A2, Ei and F 2 are given by 



(6.26) 



E 2 = -(Nf +2nN s + l)P(l+vP)-2N s 2 -2, 
F 2 = (N s + n) (n 2 + nN s + l) + n(nN s + 1) /3 2 (1 + v/3) 2 + 4N S (n 2 + nN s + l)/3(l + v/3) 



+2 (« - AT, 3 + 2N S ) p + 2N S (N s 2 + 4) 
In A 3 , E 3 and F 3 are given by 

£3 



(6.27) 



/3 2 (1 + v/3) 2 



(N s 2 + 3nN s + 3n 2 + 1)(n s 2 + 2) - 6 n 2 

+ (AN 2 + 8 + 6nN s ) (n 2 + l)/3(l + v/3) + N s 4 (3 -/3) + 6 (/3 + 3)N 2 

+12 + 3nN s + 4f3 (6.28) 

F 3 = - ((N s + nf + N s + n 3 ) (n 2 + 2) + 4 n ( 1 - n 2 ) /3 3 ( 1 + v/3) 3 
-6N s (n s 2 + niV, + 2)(n 2 + nN s + l)/3 2 (l + v/3) 2 
+ 3 08 - 3) JV. S 5 + 6 (y8 - 1) nJV. s 4 - 36JV, 3 - (24 + 15/3) nN 2 

-6(/3n 2 + 4/3 + 8)iV s - 12/8 « ,8(1 + v/3) 

/4 \ 32 

-2(3,6 + 4) -n + nN 2 - N s 5 + 8N S - 24JV. S 3 - 1(W S 5 + — n + 8nN 2 . (6.29) 

In A 4 , E 4 and F 4 are given by 

£4 = -[{2 + N S 3 (N S + 4n) + 3V, 2 (l + 2n 2 ) + 2nN s (3 + 2 n 2 ) J (N s 2 + 3) + 4nV, (l - 3n 2 ) ]/3 3 (l + v/3) 3 
— 1^36 + 2V S 3 (3V, + 8n) + 6V, 2 (5 + 2n 2 ) + 44n/V,] (n s 2 + l)/3 2 (l + v/3) 2 
+[2V, 6 08- 5) + 4niV s 5 03- 3) - 6N, 4 (/3 + 13) - 6nN s 3 (5/3 + 12) - 4/V, 2 (50 + 17/3 + 3/3n 2 ) 

-2nV, (24 + 23/3) - 36(3 + /3)j/3(l + v/3) + 4/V,, 6 08- 1) - 12/V, 4 08 + 5) - 12h/V, 3 /3 
-8V,, 2 (17/3 + 28) - 2%nN s /3 - 36(3 + 2/3) (6.30) 
[2n 4 V,, 3 + 4n 3 V,, 2 (V,, 2 + 3) + 2n 2 N s (3/V, 4 + 9N S 2 + ll) + In (n 2 + 3) (2V,, 4 + 3/V, 2 + 2) 

+/V, 7 + 6/V, 5 + 11V,, 3 + 6/V,]/3 4 (l + v/3) 4 

+8/V, (N 2 + n/Y s + l) (N 2 + n/V, + 2) (nN s + 3 + V,, 2 )/3 3 (l + v/3) 3 



F4 
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l2N s 2 (3n 3 + ll2N s 3 (4 + 3/3) + 12N S 5 (1 -j8) + SON^n 2 
+ll2N s 6 (3-0) + 6/y s 4 (26 -/3) + 30W., 2 (6 + 5/3) + 72/sjn + 4(5 -/3) A?, 7 +6(23 -/3) A', 5 
+ (310 + 46/3) W s 3 + (144/3 + 288) /Vj /3 2 (1 + r/3) 2 
8 (4 + 3M 2 ) n 2 N s /} + (8(1-3/3) W, 6 + 24 (/3 + 3) N s 4 + (280/3 + 256) W, 2 + 960)n 

+ 16(1 -/3)W S 7 + (536+ 184/3) 2V, 3 + (168 - 24/3) AT, 5 + (576/3 + 768) JV,|s(l + v/3) 
+2(l +0 2 -6p)N s 1 +(-72/3-24yS 2 + 48)W, 5 - I20n(0- 1) AT., 4 + (-8/3 2 + 168,8 + 352) 3 
+24 + jgj nW.r/3 + (768 + 6 (n 2 + 16) /S 2 + 768 p)N s + 24 f/3 + | j «y3 . 



(6.31) 



6. 2. i . ///g/i S7V7? Analysis of the Symbol Error Rate ( SER ) 

Based on (1 . 1 1), the SER of MPSK modulation can be expanded at high SNR 
using (6.21) into 

1 °° A 

Pmpsk = - V jz £ —^I M (&) (6.32) 

71 ti> (tsmpsk)^ 



where 



"r(®)= r 
Jo 



/,-(©) = I sin 2r 0d#. (6.33) 



Note that here, in contrast to (6.1 1), the exponent r is a positive integer. As such, 
(6.33) admits the following closed-form solution [50, 2.513.1]: 

As before, a first-order approximation is of key interest, giving 

/W = (G a y)~ Gd + o(r Gd ), (6.35) 
where we identify the diversity order 

G d = nN s , (6.36) 
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and the array gain 



G a = 



gMPSK 



( 



A I Gd (&)Y Gd 



(6.37) 



The result for Gd above is consistent with a result obtained via a different method 
in [8], whilst the expression for G a appears new. The high SNR results above, for 
the case /3 £ 0, are illustrated in Fig. 3. As before, the "Simulation" curves are 
based on numerically evaluating the exact SER relation (1.11), and the "Coulomb 
Fluid (Exact)" curves are based on substituting (6.2) into (1.11) and numerically 
evaluating the resulting integral. Moreover, the "Coulomb Fluid (Leading term 
only)" curves are based on (6.35), whilst the "Coulomb Fluid (Leading 5 terms)" 
curves are based on the first five terms of (6.32). Again, the leading-order approx- 
imation is shown to give a reasonably good approximation at high SNR, whilst the 
additional accuracy obtained by including a few correction terms is also evident. 

6.2.2. High SNR Analysis of the Probability Density Function ofy 

As before, based on the moment generating function expansion (6.21), apply- 
ing for fi ± 0, we can immediately take a Laplace inversion to obtain the following 
high SNR representation for the PDF of y, 



We are mainly interested in the leading order term, and so we write the above as 



Note that despite the similarity with (6.18), interestingly, these results do not co- 
incide upon taking fi — » in (6.39), due to the differences in d and A . This seems 
to indicate that the double asymptotics y — » oo and f3 — » are non-commutative. 




(6.38) 




(6.39) 
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SNR (dB) 

Figure 4: Illustration of the SER versus average received SNR (at relay) y; comparison of analysis 
and simulations. Results are shown for Nr = 2 and N s = 2, with the full-rate Alamouti OSTBC 
code (i.e., R = 1). QPSK digital modulation is assumed, such that M — 4. The relay power b is 
assumed to scale with y by setting b — \y. 



7. Characterizing A Through Painleve V 

In this section, we obtain the leading term of the large s and large y expansion 
(6.21) from a Painleve V differential equation, thus demonstrating the accuracy of 
the Coulomb Fluid approximation. We will focus in this section on the case J3 ± 0. 

Recall that the moment generating function of SNR y, regarded a function of 
s and t, given by (1.20), reads, 

/ V*\ A+ "/ J w \<i<j<n k=l \1+M + **/ 

[0,oo)' ! J 



(7.1) 



where < K n a is a normalization constant in (1.26). We also recall that c and t are 
given by 

y (1 +y)N R ~ 

t = z , where b := vy, (7.2) 



RN S 
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respectively. So as y — > oo, we note that 



Nr 
v 



A simple computation shows that M 7 (s, t) admits the following expansion for 
large ys: 



My{S, t) = 



(RN s y 



nNs 
ys 



h f n <* - ^ fl x r N '*- Xk « + 1 1 - S + • ■ • ) dXk 

•[OX,,, 1 <i<j<n k=l 



(RN s ) m ' 1 

D„[w dLag (-,o;.-^,^)](r^ D " [M;dL - ( -' u a - N " Ns)] 



R 



1 



£>«[wdLa g (-, 0, a - N s , N s )] (ys) 



nN s +l 



nD„[w dLag {;t,a- N S ,N S )] 



"■ [0i 4 !</<)<» 



(7.3) 



where WdL ag (X ?, a - A^, /V 4 .) is the deformation of the classical Laguerre weight, 
i.e., 



WdLagC*. t,a-N s ,N s ) 



e~\t + x) N \ t > 0, a - N s > -1, (7.4) 



and 'Kn a = Ait^dLagO, 0, a - N S ,N S )] is independent of 

The condition a - N s > -1 is required for the validity of the orthogonality 
relation with respect to the generalized Laguerre weight 00 (see [9]). This, 

in turn, will ensure the validity of the first two terms in (7.3), while for the multiple 
integral in the third term to converge, the condition a - N s > is required. 

For the problem at hand, a and N s are integers, satisfying a > and N s > 0. 
Therefore a > N s implies that the result for A presented below is valid for a > 0. 

Comparing the expansion (7.3) with (6.21), we see that the diversity order is 



and 



A = (RN S ) 



d = nN s , 



nNs D n [w dLag (-, t,a- N s , N s )] 
£>„[WdLag(-,0,ar- N S ,N S )]' 
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(7.5) 



(7.6) 



We see that A is up to a constant the Hankel determinant which generates a par- 
ticular Painleve V and shows up in the single-user MIMO problem studied in [13]. 
We obtain A , through a large n expansion of 

A,[>VdLag(-,f,g-A^Ag] 

A,[WdLag(-,0,(Z - N S ,N S )X 

for a > 0. We will see that A precisely matches that obtained in (6.23). 

From [13] we learned that the logarithmic derivative of D n [wd Lag (-, t, a-N s , N s )] 
with respect to t, 



H n (t) 



t— log Ai[w d Lag(-, t, a - N s , N s )] 



.-log 



A,[vVdLa g (-,f,a- A^,AQ] 
A,[w d La g (-,0,a-A^,,M)] 



(7.7) 



satisfies the Painleve V: 

>2 



K) 



[(t + 2n + a)H' n - H n + nN s f 

-4(tH' n -H n + n(n + a))(H' n )(H' n + N s ), 



(7.8) 



where ' denote derivative w.r.t. t. 

We restrict to the case where N R < N D , for which N R = n. The case where 
Nr > N D can be considered in a similar fashion as outlined in Section 6.2. Setting 
a = nj3 in the above equation, where ft = - - 1 is a fixed positive number, and 
with the change of variable 



n 

t = -, 
v 



(7.9) 



an easy computation shows that 



y„(v) := H n (n/v) 



satisfies 

,A 



(os + 2)v + i)r; + r„ - n^v,. 

+— ( - vY 'n ~ Y n + (1 +/^ 2 )(y,;)( - v 2 ^ + niV,), (7.10) 



where ' denotes derivative with respect to v. 
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We seek a solution for Y n (v) in the form 



Y n (y) = np-i(y) + p (v) + 



Pj( v ) 



7=1 



nJ 



(7.11) 



from which A (v) is found to be 

( y 

A (v) = (RN s ) nN *exp 



f 



civ 

v' 



\ CO 

/ v 



(7.12) 



(RN s ) nNs exp 



J y' ?i J v' 



(7.13) 



Substituting (7.1 1) into (7.10) leads to (7.10) taking the form 

oo 

C_2« 2 + C-\tl + C + ^ C / n~-' = 0, 
7=1 



(7.14) 



where c_ 2 depends on p-\(y) and its derivatives, and q, / = -1,0, 1,2, .. . depend 
on p-i(y), />o(v) up to p,+i(v) and their derivatives. Of course, each c, also depends 
upon y, /V. s and /3. Assuming that the coefficient of n k is zero, we find that the 
equation c_ 2 = gives us 



.((/3 + 2)y+l)// 1 (y)+/ M (y)-7V i 



4y 2 (l+/3K 1 (y)(y 2 ^i(v) + ^). 

(7.15) 



With MAPLE, the solutions of the above ODE for /?-i(y) are found to be 

( - vJ3 - 1 + ^(1 +vfif + 4v)N s 



P-i(v) 
P-i(v) 



2v 



(-VJ3-1- +vfif + 4v)N s 



2v 



(7.16) 
(7.17) 



P-i(v) = (2 + J3 + l - j C 2 + N s + 2 V(l + /« + C 2 ), (7.18) 
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where C 2 is a constant of integration 7 . 

The equation c-\ = is a coupled differential equation involving both p~\{v) 
and po(v) given by 

2v 3 (Ki(v)) 2 + 2vY_ 1 (v)/m(v) + v(l - vyS) // (v) + po(v)| (n s - v/»_i(v)) 



-V 



(/?V + 2 OS + 2) v 2 + v y (v) - 2 vVi(v) 2 + (l + OS + 2) v)p (v) 



-2 v 5 (// jCv)) 3 = (v (v + 1) OSv + 1) p|,(v) - (v - 1) po(v))p_i(v). 



(7.19) 



With p-i(v) given by (7.16), chosen to match the result from the Coulomb Fluid 
(6.23), we find that the first order ODE in po(v) has the solution, 

(2 + p + v/? 2 - p J (I + vpf + 4vW 

Po(v) = - 7 — 2 r . (7.20) 

2((1 + vpf +4v) 

We disregard the second and third solutions for p~i(y); as these would lead to po(y) 
which do not generate the A in agreement with that obtained from the Coulomb 
Fluid method. 

Substituting p-\(y) from (7.16) and po(v) from (7.20) into c = gives, 
N s v 2 [(p 2 (2 +0) v 2 + 2 (j3 2 + 2p + 2) v + 2 + f3)N 2 - (1 +/?) v] 



pi(y) 



((1 + vpf + 4v) 



5/2 



A^V(1 + (2+ P)v)p 



((1 +vyS) 2 + 4v)~ 
Hence, A (v) has a large n expansion, 



(7.21) 



A (v) 



?o(v) 



1 + 



<7iO) 



+ 01- 



where 



<7o(v) 



(1 + v/?) 2 + 4v) 2 ' " (l + 2v + v/J + J(l + v/?) 2 + 4v)~ 



v ^( (1+Vy8) 2 +4v )^(l +y 6) ¥(2+ V 
exp (-^(l + v/? - V(l+v/?) 2 + 4v)) . 



, (RNs) 

' 2^(2n+iV») 



(7.22) 



M2+/S) 



(7.23) 



7 Note that the constants of integration in (7.16) and (7.17) are zero. 
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The leading term of A in (7.22) agrees precisely with the A computed via the 
Coulomb Fluid method in (6.23). 

Using a method similar to the cumulant analysis of Section 5, we can also 
compute the first correction term to A (i.e., the quantity qi(v)) using (7.21), which 
reads 



qi(y) = 



N s ( (2N S 2 - l)y6 2 + 2 (M s 2 - l) (1 + £>) 



24(1 + /?) 



1 _ v^£v + 3/3 + 6) 
^"((l+y^) 2 +4v) 3/2 



(2/3 v + 4 v + 1) N s 3 (2 + )8) (2W, 2 - l) N s 



2/i((l + v/3) 2 + 4v) 8 ((1 + y/3) 2 + 4y) 3/2 (1 + fi) 
(2 + /3)A^y 2 



(2+/3)y+- 



\3/2 ' 



((1 +y/3) 2 +4y) 

Higher order corrections could also be obtained in a similar way. 



(7.24) 



8. Conclusion 

In this paper, we have introduced two methods for characterizing the received 
SNR distribution in a certain MIMO communication system adopting AF relay- 
ing. We showed that the mathematical problem of interest pertains to computing a 
certain Hankel determinant generated by a particular two-time deformation of the 
classical Laguerre weight. By employing the ladder operator approach, together 
with Toda-type evolution equations in the time variables, we established an exact 
representation of the Hankel determinant in terms of a double-time PDE, which 
reduces to a Painleve V in various limits. This result yields an exact and funda- 
mental characterization of the SNR distribution, through its moment generating 
function. Complementary to the exact representation, we also introduced the lin- 
ear statistics Coulomb Fluid approach as an efficient way to compute very quickly 
the asymptotic properties of the moment generating function for sufficiently large 
dimensions (i.e., for sufficiently large numbers of antennas). These results, which 
have only started to be used recently in problems related to wireless communica- 
tions and information theory, produced simple closed-form approximations for the 
moment generating function. These were employed to yield simple closed-form 
approximations for the error probability (for a class of M-PSK digital modula- 
tion), which were shown via simulations to be remarkably accurate, even for very 
small dimensions. 
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To further demonstrate the utility of our methodology, we employed our asymp- 
totic Coulomb Fluid characterization in conjunction with the PDE representation 
to provide a rigorous study of the cumulants of the SNR distribution. Starting with 
a large-/? framework, we computed in closed-form the finite-^ corrections to the 
first few cumulants. It was seen that the Coulomb Fluid approach supplies the cru- 
cial initial conditions which are instrumental in obtaining asymptotic expansions 
from the PDE. We also derived asymptotic properties of the moment generating 
function when the average SNR was sufficiently high, and in such regime ex- 
tracted key performance quantities of engineering interest, namely, the array gain 
and diversity order. 
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Appendix A. Characterization of Hankel Determinant Using the Theory of 
Orthogonal Polynomials 

In this section, we describe the process by which we characterize the Hankel 
determinant by using the theory of orthogonal polynomials and their associated 
ladder operators. See [9, 13, 14, 16, 19] for the background to this theory. 

Consider a sequence of polynomials {P n (x)} orthogonal with respect to the 
weight function w AF (x, T, t), given by (1.21) 



w AF (x, T, t) = 




I , 0<x<oo, (A.l) 



i.e., 

oo 

J P n (x)P m (x)w AF (x, T, t)dx = h n S n , m , (A.2) 
o 

where h n is the square of the L 2 norm of P„(x). The Hankel determinant (1.1) is 
reduced to the following product: 

n-l 

DAwaf] = Y]hj. (A3) 

7=0 

Our convention is to write P n (x) as 

P n (x) := x" + PjOi)*"" 1 + p 2 (n)x n ~ 2 + ■■■ + Pj(0). (A.4) 

Hence, this implies that properties of the Hankel determinant may be obtained 
by characterizing the class of polynomials which are orthogonal with respect to 
w AF (x, T, t), over [0, oo). It is clear that the coefficients of the polynomial P n (x), 
Pi(n), will depend on 7, J, a and N s ; for brevity, we do not display this dependence. 
From the orthogonality relation, the three term recurrence relation follows: 

xP n (x) = P n+1 (x) + a n P n (x) + p n P n -x{x), n = 0, 1, 2, . . . (A.5) 

with initial conditions 

P (x) si, and 0oP-i(x) = 0. 

The main aim is to determine these unknown recurrence coefficients a n and J3 n 
from the given weight. Substituting (A.4) into the three term recurrence relation 
results in 

a n = Pi(ra) - Pi(« + 1) (A.6) 
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where p^O) := 0. Taking a telescopic sum gives 



n-1 

Y^j = -PiOO- (A.7) 

;=o 

Moreover, combining the orthogonality relationship (A. 2) with the three term re- 
currence relation leads to 

Pn = 7^, (A.8) 

h n -\ 

which can also be expressed in terms of the Hankel determinant D n in (A.3) 
through 

fin = ^± (A.9) 

since 

, D n+1 

n„ = . 

D 

Ladder Operators, Compatibility Conditions, and Difference Equations 

In the theory of Hermitian random matrices, orthogonal polynomials plays 
an important role, since the fundamental object, namely, Hankel determinants or 
partition functions, are expressed in terms of the associated L 2 norm, as indicated 
for example in (A.3). Moreover, as indicated above, the Hankel determinants 
are intimately related to the recurrence coefficients a n and fi n of the orthogonal 
polynomials (for other recent examples, see [19, 51-53]). 

As we now show, there is a recursive algorithm that facilitates the determi- 
nation of the recurrence coefficients a n and f3 n . This is implemented through the 
use of so-called "ladder operators" as well as their associated compatibility condi- 
tions. This approach can be traced back to Laguerre and [54]. Recently, Magnus 
[11] applied ladder operators to non-classical orthogonal polynomials associated 
with random matrix theory and the derivation of Painleve equations, while [10] 
used the associated compatibility conditions in the study of finite n matrix mod- 
els. See [14, 16, 19] for other examples of the application of this approach. 

From the weight function w AF (x, T, t), one constructs the associated potential 
v(;c) through 

v(jc) = -logw AF (jc) = x-alogx-N s logl^-^-), (A.10) 

\T + x) 
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and therefore, 



a 



v'(x) = 1 - - 

x x + t x + T 



+ 



(A.11) 



As shown in [15], a pair of ladder operators, to be satisfied by our orthogonal 
polynomials of interest, are expressed in terms of v(x) and are given by 



d 

— + B„(x) 
ax 



— - B n (x) - v'(x) 
ax 



P n (x) =B n A n (x)P n -i(x), 

Pn-l(x) = ~ A n ^( X )P n (x), 



(A.12) 



where 



CO 

K J x-y 



B n (x) 



h 



i r v'(x) - v'O) 

— P n (y)Pn-i(y)w A F(y)dy, 

n-i J x-y 



(A.13) 



and where, for the sake of brevity, we have dropped the t, T dependence in waf- 

Moreover, there are associated fundamental compatibility conditions to be sat- 
isfied by A n (x) and B n {x), which are given by [16] 

B n+1 (x) + B n (x) = (x- a n )A n (x) - v'(x), (S i) 

1 + (x - a n )[B n+l (x) - B„(x)] = B n+l A n+l (x) -8 n A n -i(x). (S 2 ) 

These were initially derived for any polynomial v(x) (see [55-57]), and then were 
shown to hold for all x 6 C U {oo} in greater generality [16]. 

We now combine (S i) and (S 2 ) as follows. First, multiplying (5 2 ) by A n {x), it 
can be seen that the RHS is a first order difference, while (x-a„)A n (x) on the LHS 
can be replaced by B n+l (x) + B n (x) + v'(x) from (Si). Then, taking a telescopic 
sum with initial conditions 



B (x) = A_j(x) = 



leads to the useful identity 



YjAjix) + B l ?(x) + v\x)B n (x) = 8 n A n (x)A n ^(x). 
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The condition (S' 2 ) is of considerable interest, since it is intimately related to the 
logarithm of the Hankel determinant. In order to gain further information about 
the determinant, we need to find a way to reduce the sum to fixed number of 
quantities; for which, (S' 2 ) ultimately provides a way of going forward. 

Remark 6. Since our V(x) is a rational function of x, we see that 



is also a rational function of x, which in turn implies that A n (x) and B n (x) are 
rational functions of x. Consequently, equating the residues of the simple and 
double pole at x = 0, x = -T, x = -t on both sides of the compatibility conditions 
(Si), (S2) and (S' 2 ), we obtain equations containing numerous n, T and t depen- 
dant quantities; which we call the "auxiliary variables " ( to be introduced below). 
The resulting non-linear discrete equations are likely very complicated, but the 
main idea is to express the recurrence coefficients a n and B n in terms of these 
auxiliary variables, and eventually take advantage of the product representation 
(A3) to obtain an equation satisfied by the logarithmic derivative of the Hankel 
determinant. 

Now substituting (A. 14) into (A. 13) which define A n (x) and B n (x), and fol- 
lowed by integration by parts, we find 



V(x) - V(y) 




(A. 14) 



x-y 



xy (x + t)(y + t) (x + T)(y + T)' 



A„(x) = 



B n (x) = 




(A.15) 



(A.16) 



where 



00 



00 













00 



00 













(A. 17) 



are the auxiliary variables. 
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Difference Equations from Compatibility Conditions 

Inserting A n (x) and B n {x) into the compatibility conditions (S i), (S 2 ) and (S' 2 ), 
and equating the residue as described above yields a system of 12 equations. Note 
that there are actually 14 equations, the compatibility conditions (Si) and (S 2 ) 
also have 0(x°) terms which yield equations that lead to = 0. 

The compatibility conditions (Si) give the following set of equations, 

a - a n (R n + 1 - K), (A. 18) 
N s -(a n + t)R* n , (A.19) 
-N s + (a n + T)R n . (A.20) 

Similarly, the compatibility condition (S 2) gives the following set of equations, 

- a n (r n +i - r„ - 1 - r* +1 + f n ) = jff„+i ~ /?« 

+J3 n+1 (R n+1 - R* n+1 ) - p n (R n -x -K-i)XA.21) 
a + a„)(r;-r; !+1 ) = p n+x R* n+l -p n R* n _ x , (A.22) 
(T + a„)(r„ - r„ +1 ) = fi n +iR n +i -fi n R n -i- (A.23) 



^+1 + r n - r n+l -r n -2n-\ = 



To obtain a system of difference equations from the compatibility condition (S' 2 ) is 
somewhat more complicated, but carrying out the process, equating all respective 
residues in (5^) yields six equations. The first three are obtained by equating the 
residues of the double pole at x = 0, x = -t and x = -T respectively: 

- 2r n r* - (2n - N s + a)r n 
+(2n +N S + a)r* n + n(n + a) = /3 n - f$ n (R n R* n _ x + 

+MRn + Rn-l) ~ Pn(K + R^J, (A.24) 

r* n K-N s ) = /3 n R* n RU, (A.25) 
r n (r n -N s ) = p n R n R n -u (A.26) 

while the last three are given by equating the residues of the simple pole at x = 0, 
x = -t and x = -T respectively: 

Vrs D n , (2r:-N s )(r n -n-r:)-ar: r (N s - 2r„)(r„ - n - Q + ar„ 

2JRj - Rj) + + + r„ - r n 

7=0 



'^Mflrp-B . P P > ^ PniK+K-i- 2 KK-d f3„(R„+R„- l +2R„R„- l ) 
- + -\p n {R n R n -\ +R„R n _ l ) + , 



(A.27) 
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(1-1 



I*;- 

7=0 



(2r* n - N s )(r„ -n- r* n ) - ar* (N s - 2r n )r* n + N s r n , 
1 1- r„ 



T-t 



T-t 

Pn(R* n Rn-l + R n R* n -l)> (A.28) 



\n n - 2r„)(r„ - n - r*) + ar n (N s - 2r n )r* n + N s r n 

> Rj + + +r„ 

L-i 3 T T-t 

7=0 

= j n{ Rn+R n -^2R n R n -i) + ^ _ J_y (KRn ^ + ^ i} . (A . 29) 

Remark 7. can that (A.27) is a combination of (A.28) and (A.29), and 

serves as a consistency check. We also see that (A.28) and (A.29) respectively 
gives us the value ofYijR) and Y,jRj automatically in closed form. Hence we 
can also obtain any linear combination of these sums in closed form, which is a 
crucial step in obtaining a link between a linear combination of the logarithmic 
partial derivatives of the Hankel determinant, and the quantities /3,„ r n and r*. 
This step would not be possible without (S ' 2 ), also known as the bilinear identity. 

Analysis of Non-linear System 

Whilst some of the above difference equations (A. 1 8)-(A.29) look rather com- 
plicated, our aim is to manipulate these equations in such a way as to give us 
insight into the recurrence coefficients a n and /3„. Our dominant strategy is to al- 
ways try to describe the recurrence coefficients a n and /3„ in terms of the auxiliary 
variables R n , r n , R* n and r*. 

The first thing we can do is to sum (A. 18) through (A. 20) to get a simple 
expression for the recurrence coefficient a„ in terms of the auxiliary variables R„ 
and R* n : 

a„ = TR n - tR* n + a + 2n + 1 . (A.30) 

We can immediately see that by taking a telescopic sum of the above from j = 
to j = n - 1, and recalling equation (A.7), gives rise to 

n-l n-l n-l 

TY J R i- t Yu r i +n{n+a) = Z a > 

7=0 7=0 7=0 

= - Pl (n). (A.31) 



60 



Looking at (S 2 ) next, we sum equations (A.21) through (A. 23) to get 

a n = fin+i ~fin + T(r n - r n+ i) + t(r* n+1 - r* n ). (A.32) 

Taking a telescopic sum of the above from j = to j = n - 1 and rearranging 
yields 

p n = Tr n - tr* n - Pl (n). (A.33) 

We are now in a position to derive the following important lemma, which de- 
scribes the recurrence coefficients a n and /?„ in terms of the set of auxiliary vari- 
ables: 

Lemma 2. The quantities a n and/3 n are given in terms of the auxiliary variables 
R n , r n , R* n and r* as 

a„ = TR n - tR* n + a + In + 1 , (A.34) 

Pn = (1 + Rn) ^ (Rn - 1) 2r„r„ 

1 + K„ - K n I K n K n 

-(In - N s + a)r n + (In + N s + a)r* n + n(n + a) . 

(A.35) 

Proof. Equation (A.34) is a restatement of equation (A. 30). 

To obtain (A.35), we use the equations obtained from^). We eliminate 
and i?„_! from (A. 24) using (A. 25) and (A. 26) respectively, and then rearrange to 
express fi„ in terms of the auxiliary quantities. □ 

Toda Evolution 

In this section, n is kept fixed while we vary two parameters in the weight 
function (1.21), namely T and t. The other parameters, a and A^ 4 are kept fixed. 
Differentiating (A. 2), w.r.t. T and t, for m = n, gives 

d T (\ogh n ) = -R n , (A.36) 
d t (\ogh n ) = R* n (A.37) 

respectively. Then, from equation (A. 8), i.e. /3„ = h n jh n -\, it follows that 

d T /3 n = fi„(Rn-l ~ Rn), (A.38) 

dfi n = MK-K-il (A.39) 
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Applying d T and d t to the orthogonality relation 

CO 

J w AF (x, T, t)P n (x)P n -i(x)dx = 



results in the following two relations: 

drPM) = r n , (A.40) 
d, Pl (n) = -r* n . (A.41) 

Note that in the above computations with d t and d T we must keep in mind that 
the coefficients of P n (x) depend on t and T. 
Using the above two equations, we get that 

(Td T + tdMn) = Tr n -tr* n , 

= p l (n)+fi n , (A.42) 

where the second equality follows from (A. 33). 

Using equation (A. 6), i.e. a n = p^rc) - Pj(n + 1), we get that 

d T a n = r n -r n+y , (A.43) 
d t a n = r: +l -r* n . (A.44) 

Now, combining (A.43), (A.44), (A.33), and (A.6), and similarly using (A.38) 
and (A. 39), we arrive at the following lemma: 

Lemma 3. The recurrence coefficients a n andfi n satisfy the following partial dif- 
ferential relations: 

(Td T + td t - l)a n = fi n -Pn+u 
(Td T + td t ~ 2)yS„ = M^n-i ~ & n ). 

In the second of the above two equations, we eliminate a n using the first equation 
to derive a second order differential-difference relation for f3 n : 

(T 2 d 2 TT + 2Ttd 2 Tt + t 2 dl) log A, = Pn-\ ~ 2/3 n +J3, 1+l - 2, 

which is a two-variable generalization of the Toda molecule equations [58]. 
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Toda Evolution ofHankel Determinant 

Recall that the moment generating function is related to the Hankel determi- 
nant through equation (1.23): 

1 / T \ n ^ s 



l tl V 
My(T, t) = - - D n (T, t), 



while the Hankel determinant is related to h„(T, t) by the relation D n (T, t) = 
W%lhj{T,t). 

We compute the partial derivatives of the Hankel determinant by taking a tele- 
scopic sum from j = to j = n-l of the partial derivatives of log h n , i.e. equations 
(A.36) and (A.37). Using (A.3), we then obtain 

n-l 

d T (\ogD n ) = -J]Rj, (A.45) 

7=0 
n-l 

d t (logD„) = J]R*. (A.46) 

7=0 

Hence we now have expressions of the partial derivatives of D n (T, t) in terms of 
the auxiliary variables R n and R* n . 
At this point, recall that 

H n (T,t) = (Td T + td t )\ogD n (T,t). (2.1) 

This quantity related to the logarithmic derivative of the Hankel determinant is 
of special interest as our goal is to find the PDE that H„(T, t) satisfies. From the 
definition of H,,(T, t), along with (A.45) and (A.46), we find that 

n— 1 n~\ 



H n (T,t) = -Tj^Rj + tJ^R*!, (A.47) 

7=0 7=0 

= p t (n) + n(n + a), (A.48) 



where the second equality is a result of (A.31). 

At this point, we may express the fundamental quantity H n (T, t) in terms of 
the auxiliary variables and /3„. This is given in the following key lemma. 



63 



Lemma 4. 

H n (T, t) = 2r n r* n + (2n - N s + a + T)r„ - {In + N s + a + t)r* 

_(1 + R n f» {r "- Ns) + (1 - K) rnirn ~ Ns) +/3 n R n -f3 n R* n . 

K n K n 

(A.49) 

Proof. Using (A.28) and (A. 29), we obtain the sum 

n— 1 n— 1 

in closed form. All that remains is to eliminate R* n l and R n -\ from this equation. 
We do this using (A.25) and (A.26). □ 

To continue, we apply td T + td, to equation (A.48), and find that 

(Td T + td t )H n = pM+Pn, (A.50) 

where we have used (A. 42) to replace (Td T + td t )p x (ri) by Pj(n) + yS„. 

We see that equations (A.48) and (A.50) can be regarded as a set of simulta- 
neous equations for p { {ri) and /?„. Solving for p^n) and /?,„ we find that 

Pj(n) = H n - n(n + a), 

and 

p n = (Td T + td,)H n - H n + n(n + a). (A.51) 

Before completing the proof of Theorem 1, we show here that our Hankel 
determinant D n (T, t) is related to the r-function of a Toda PDE. 

Substituting the definition (2.1) into (A.51) and together with (A. 9), we find 



{Td T + td t f log A, - (Td T + td t ) \ogD n + n(n + a) = D " +i ^ nl 

D n 



Defining 



n(n+a) 



D n (T,t):=(Tty—D n (T,t), 
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after some computations, D n (T, t) satisfies the following equation: 



{^d\ T + 2d\ t + L d A log£ ) n {T,t) = (A.52) 

n 

The above equation is a two parameter generalization of the Toda molecule equa- 
tion [58], and hence we identify D„(T,t) as the corresponding r-function of the 
two-parameter Toda equations. 

A reduction may be obtained by writing D n (T, t) as 

n(n+a) 

>T\— 



D„(r,f) = (y) 2 ® n (T), 



and equation (A.52) is reduced to a 1 -parameter Toda equation, 



^ r io g o„(r) 



Partial Differential Equation for H n (T, t) 

We differentiate (A. 48) with respect to T and t, and make use of the expres- 
sions for drp^n) and d t p y {ri), i.e., (A.40) and (A.41) respectively, to give 

d T H n = r n , (A.53) 
d,H n = -r*. (A.54) 

Thus, we now have expressed r n and r* in terms of the partial derivatives of H n . 

Next we derive representations of R n and R* n in terms of H n and its partial 
derivatives. 

The idea, following from previous work [19, 51, 52, 59], is to express R„ and 
1 /R n in terms H n and its derivatives w.r.t. t and T. Similarly, for i?* and 1 /R* r . 
We start by re-writing (A. 26) as 

a D r n (r n - N s ) 

Pn -fVi-1 - 



R 



and substitute into (A.38) to arrive at 



a a r n( r n N s ) 

9tPu = ~ PnRn, 



which is a linear equation in R n and 1 /R n . 
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Substituting r n given by (A. 53) and J3 n given by (A. 51) into the above equation 
produces 

K n 

-{T(d T H n ) + t(d,H n ) - H n + n(n + a))R n . (A.55) 
Going through a similar process we find, using (A. 25), (A.51) and (A. 54), 

. _ (W(^> 

+(T(d T H n ) + t(d t H n ) - H n + n(n + or))^. (A.56) 

The above two equations are quadratic equations in R n and R* n . Solving for them 
leads to 

2(T(d T H n ) + t(d t H n )-H n + n(n + aj)R n = -(Td 2 TT + td 2 Tt )H n ± A,{H n ), 

(A.57) 

2(T(d T H n ) + t(d t H n ) -H n + n(n + aj)R* n = (Td 2 Tt + td 2 tt )H n ± A 2 (H n ), 

(A.58) 

where Ai(H n ) and A 2 (H„) are defined by (2.4) and (2.5) respectively, where we 
have left out the notation that indicates that Ai and A 2 also depends upon the 
partial derivatives of H n . 

In the last step we substitute (A.53), (A.54), (A.51), R n given by (A.57) and R* n 
given by (A.58) into (A. 49). After some simplification, we obtain the PDE (2.3), 
completing the proof of Theorem 1 . 

Appendix B. Some Relevant Integral Identities 

The following integral identities, taken from [13], are referenced for use in the 
Coulomb fluid derivations. 
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Let W(t) = VO + a)(t + b): 

b 



I 



dx 



I 



yl(b - x)(x - a) 

a 

dx 

(x + t) y/(b - x)(x - a) 



= n, 



D 



xdx 



V(£> - x)(x - a) 



y/(b-y)(y-a) 



p _^ ^ , 

(y-x)(y + t) 



D 

2^1 



y/(b — x)(x - a) 

log(* + 
x V(& - x)(x - a) 



2^1 



y/(b - x)(x - a) 



dx 



a 

2^1 



log(x + t) 



(x + t) y/(b - x)(x - a) 



zdx 



n 



w(ty 

(a + b) 



x + t 



i r iogu + o „. _ 

2n J yl{b - xYx - a) 



log 



' (y/ab + W(t)f - t 2 " 

2 



(B.l) 



(B.2) 



(B.3) 



(B.4) 



(B.5) 



, (B.6) 



2 

jdogo + o , ( VF+^ - VF+75)" 



, (a + b) 
+ — - — log 



'(W{t) + if -ab 



4t 



(B.l) 



— log! — -L= + — ^=|,(B.8) 



f — 

J x^a' 



dx 



x 2 + b'x + c' 



1 , ( 2c' + b'x + 2 Vc 7 Va'x 2 + + c' 
: log 



(B.9) 



Vc 7 

where (c' > 0), 

logf V where (c > 0,b' 2 = Adc). (B.10) 

yTj' \2c' + fc'x/ 
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Finally, note the well-known parameterization 



log(A + 5) = logA + 



Bdri 
A + rjB 



(B.ll) 



Lemma 5. From the above integrals, we have 



- [— 

2nJ (jc + 



log(x + 1 2 ) 



t\) yj(b - x)(x - a) 



-_dx 



1 



2W(h) 



log 



\w{ty) + W{t 2 )f -{h -t 2 f 

2 

( yjt\ + a + yftT+b) 



(B.12) 



Proof. We rewrite log(x + t 2 ) using the Schwinger parameterization (B.ll) and 
then use (B.2) to get 



1 b 



LHS 



(5.12) 



log h + r r 

+ a)(ti + b) 2tt J J 

o a 



xdxdn 



(t 2 + x?])(x + t\) - x)(x - a) 



Now we use the partial fraction decomposition 

x t\ 



if 2 + xrj)(x + h) (hrj- t 2 ){x + h) (hT]- t 2 ){t 2 + xt]) 
and then integrate over the x variable using (B.2) to get 

\ogt 2 



LHS 



(5.12) 



2 V(?i + «X*i + b) 
i 



2 J ?H7- 



*2 



?2 



log(f2-fl) 



V(*i + flX^i + 5) y(? 2 + «77)(f 2 + fej) J ' 



lie 



2 V(fi + + 5) 



-Ilhn f 

2^0 J ( flJ? _ 



(B.13) 



In our problem ? 2 = *i = 7" or f 2 = 7", h = f, hence f 2 /fi = (f /T') ±l = 
(t/T) ±l = (1 + cs) ±l . For the case s = 0, there exists another pole within the 



lie 



integrand, and so we replace the j . . . with f ■ ■ ■ so that we may invoke (B.9). 
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To evaluate the remaining integral in (B.13), we first make the change of vari- 
able y = t\t] - t 2 , giving 



lum f ^ = iHm f ^ 

2*-+oJ (til] - t 2 ) J(t 2 + arf)(t 2 + brf) 2 <^o J y J^jfa + a) + ay)(f 2 (f 1 + + fey) 

— ?2 

(B.14) 

Within the square root term, we have a quadratic function in y, given by 

aby 2 + ((a + b)t\t 2 + 2abt 2 )y + ffOi + a)0i + b). 

For our problem, where a and £ are given by (3.25) and /3 > 0, we see that the 
discriminant of the quadratic form is given by 

t 2 (b - af > 0, 

while the constant term is 

+ 2ti(2 + JZ) + p 2 ) >0. 
Hence we may invoke (B.9) and then take the limit e — » to get 



/ 2 V(fi + + b) V(f 2 + + b) + 2ab + (t 2 + h)(a + b) + 2t 2 h \ \ 

+ log = . (B.15) 

\ 2 V(?i + + + 2ti + a + b II 

Substituting this back into (B.13) and after some algebra, we get (B.12). □ 

Appendix C. Differential Equations For Large Scale Corrections To Cumu- 
lants 

Appendix C.l. K 2 (t') 

The correction terms fi{?) for i = 1,2,3,4 satisfy the following differential 
equations: 



((;' + p) 2 +4t'f I df l (f) 
2?{\+p) df 
((f + p) 2 +4t'fdf 2 {t') 
2N s f{\ + p) df 

((f +P) 2 + 4t') l5/2 dMt') 
2f{\+p) df 

{{f + pf+4t'f df A {f) 
2N s f(l+P) df 



£V), (C.l) 

4 2 V), (C2) 

l 2 \t'), (C.3) 

lf{t'). (C.4) 
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The functions ^ (f), i = 1, 2, 3,4 are given by 

4V) - 3?' 4 - 3 + 2) f' 3 - 2 (6y8 2 +y? + l) f' 2 - 3 + 2) /3 2 ?' + 3/3 4 . (C.5) 

l (2 \t') = -I6t' 6 -2 (J3 + 2) t' 5 +3(23 p 2 -$J3-$)t' 4 

+4(J3 + 2)(l5/3 2 - 2/3 - 2) t' 3 - 2j3 2 (lj3 2 - 48/3 - 48) t' 2 

-\Sj3 4 (j3 + 2)t' + tf>. (C.6) 

lf{t') = 80f' 8 - 60 (/3 + 2) ?' 7 - 3 (217yS 2 - 47 yS - 47) f' 6 

-08 + 2) (497yS 2 -55/3-55) ?' 5 

+3(l75/3 4 -499 /3 3 - 481/3 2 + 36/3 + 18)?' 4 

+ 15(yS + 2) (42/3 2 - 29/3 - 29) f' 3 /3 2 

+5/3 4 (7/3 2 + 192/3 + 192)f' 2 -81(/3 + 2)/3V + 3/3 8 . (C.7) 
4 4) (f') - -540f' 10 + 304(/3 + 2)f' 9 + 4(l437/3 2 -496/3-496)f' 8 
+36(/3 + 2) (185/3 2 - 44/3 - 44) t' 1 
- (4593^ 4 - 30152 /3 3 - 26584/S 2 + 7136/3 + 3568) ?' 6 
-4(/3 + 2) (2751^ - 4658/3 3 - 4442/3 2 + 432/3 + 216) t' 5 
-3J3 2 (1227 Z3 4 + 10456^ + 3944 /3 2 - 13024/3 - 6512) ?' 4 
+36/3 4 08 + 2) (53y3 2 - 302/3 - 302) t' 3 

+3/3 6 (275/3 2 + 1632/3 + 1632) ?' 2 - 108y6 8 (fi + 2) t' + /3 10 . (C.8) 



Appendix C.2. Kj,(t') 

The correction terms for i = 1, 2, 3, 4 satisfy the differential equations 



((;'+y3) 2 +4Q 11/2 Jg 1 (Q 
6f(l+/3) A' 
((f'+/3) 2 +4Q 7 </ g2 (f') 
6N s t'{\+p) df 
((t'+P) 2 + 4t') 17/2 dgi(t') 
6f(l+/3) 
((f / +^) 2 + 4O 10 Jg 4 (f0 
6/V 4 ?'(l + /3) <i?' 



i?V), (C.9) 

lf{t'), (CIO) 

4 3 V), (en) 

#V). (C.12) 
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The terms lf(f), i = 1, 2, 3, 4 are given by 

4V) = -t' 6 +9(j3 + 2)t' 5 +3(5/3 2 -6j3-6)t' 4 -(\5/3 2 + 2j3 + 2)(j3 + 2)t' 3 

-2 (15 yS 2 + 4yS + 4) ?' 2 /3 2 -608 + 2) ^ 4 + 4/3 6 + ?'/V s 2 [ - ?' 5 + 08 + 2) ?' 4 

+7f' 3 y 6 2 + 5 (/i + 2) f' 2 /3 2 - 2 (y3 2 - 4,6 - 4) t'/3 2 - 2 (fi + 2) ,6 4 ] . (C. 13) 

4 2 V) - 16?' 8 -61(/3 + 2)f' 7 - 8 (28/3 2 - 3/3 - 3)f' 6 -7(/3 + 2)(3/3 2 +4/3 + 4)?' 5 

+2 ( 189/3 4 - 76/3 3 - 92 /3 2 - 32/3 - 16) t' 4 + lp 2 (p + 2) (39 p 2 - 4/3 - 4) t' 3 

-4J3 4 (7J3 2 - 90 p - 90) f' 2 - 27^ 08 + 2) ?' + 2/3 8 . (C.14) 
ifV) - -80f' 10 + 320(/3 + 2)?' 9 + 5(253/3 2 - 159/3- 159)?' 8 

-06 + 2)(311/8 2 - 125/3- 125)?' 7 

-2 ( 1946 Z3 4 - 1015/3 3 - 1062/? 2 - 94/3 - 47) t' 6 

-2 06 + 2) ( 1498 Z3 4 - 413/3 3 - 422 j3 2 - 18,3 - 9) t' 5 

+p 2 (1085 Z3 4 - 7567 /3 3 - 6843 /3 2 + 1448 p + 724) f' 4 

+5P 4 08 + 2) (341 /3 2 - 395/3 - 395) f' 3 

+2/3 6 (91 p 2 + 1230/3 + 1230) t' 2 - 142 /3 8 (fi + 2)T + 4/3 10 

+?'/V. v 2 [ - 50f' 9 + 70 08 + 2) f' 8 + 2 (325^ 2 - 97/3 - 97) ?' 7 

+2 06 + 2) (305 p 2 - 53 p - 53) t' 6 

-2 (380/3 4 - 1535/3 3 - 1407/? 2 + 256,3 + 128) t' 5 

-2 06 + 2) (700y3 4 - 981 fi 3 - 949 p 2 + 64p + 32) f' 4 

-2^ 2 (215/3 4 + 1982/3 3 + 918/3 2 - 2128/3 - 1064) ?' 3 

+2P 4 08 + 2) (125 /3 2 - 684,6 - 684) t' 2 

+22/3 6 (5/3 2 + 28/3 + 28) f' - 10/3 8 (/3 + 2) ]. (C. 15) 

4 4 V) - 1 080 f' 12 - 3460 03 + 2) f' 11 - 32(619/3 2 - 338/3 - 338) ?' 10 
-2 06 + 2) (2045 p 2 + 2Q4P + 204) t' 9 
+20(3309/3 4 - 4120/S 3 - 4040/3 2 + 160/3 + 80) ?' 8 
+ 03 + 2) (81555/3 4 - 62404y3 3 - 62084y3 2 + 640/3 + 320) t' 1 
-80/3 2 (57/3 4 - 3845^ 3 - 2401 p 2 + 2888/3 + 1444) t' 6 
-3p 2 06 + 2) (l 8777 Z3 4 - 58336/3 3 - 50528/3 2 + 15616^ + 7808) f' 5 
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-6/J 4 (3915/? 4 + 24728/3 3 - 2792 p 2 - 55040,3 - 27520) t' 4 
+45/3 6 (/i + 2)(l09/3 2 - 1172/3 - 1172) f' 3 

+4yS 8 (725 /3 2 + 4104,8 + 4104) t' 2 - 275p 10 (J1 + 2)t' + 2p 12 . (C.16) 
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